Perfectly Matched Layers and High Order Difference Methods for Wave Equations

Dedicated to my father

Mr. Ambrose Duru (1939 – 2001)

List of papers

This thesis is based on the following papers, which are referred to in the text by their Roman numerals.

I K. Duru and G. Kreiss, (2012). A Well–posed and discretely stable perfectly matched layer for elastic wave equations in second order formulation, Commun. Comput. Phys., 11, 1643–1672

(DOI:10.4208/cicp.120210.240511a).

contributions: The author of this thesis initiated this project and performed all numerical experiments. The manuscript was prepared in close

collaboration between the authors

II G. Kreiss and K. Duru, (2012). Discrete stability of perfectly matched layers for anisotropic wave equations in first and second order formulation, (Submitted).

contributions: The author of this thesis initiated this project and performed all numerical experiments. The manuscript was prepared in close

collaboration between the authors.

III K. Duru and G. Kreiss, (2012). On the accuracy and stability of the perfectly matched layers in transient waveguides, J. of Sc. Comput., DOI: 10.1007/s10915-012-9594-7. In press.

contributions: The author of this thesis initiated this project performed all numerical experiments and had the responsibility of writing the paper. The remaining time was spent between the author and his advisor correcting misconceptions, improving the texts and the theory.

IV K. Duru and G. Kreiss, (2012). Boundary waves and stability of the perfectly matched layer. Technical report 2012-007, Department of Information Technology, Uppsala University, (Submitted)

contributions: The author of this thesis initiated this project and had the responsibility of writing the paper. The remaining time was spent between the author and his advisor correcting misconceptions, improving the texts and the theory.

V K. Duru and G. Kreiss, (2012). Numerical interaction of boundary

waves with perfectly matched layers in elastic waveguides. Technical

report 2012-008, Department of Information Technology, Uppsala

University, (Submitted).

contributions: The author of this thesis initiated this project, performed all numerical experiments and had the responsibility of writing the paper. The remaining time was spent between the author and his advisors correcting misconceptions, improving the texts and the theory.

VI K. Duru, G. Kreiss and K. Mattsson, (2012). Accurate and stable boundary treatments for elastic wave equations in second order formulation, (Submitted).

contributions: The author of this thesis initiated this project, performed all numerical experiments and had the responsibility of writing the paper. The remaining time was spent between the author and his advisors correcting misconceptions, improving the texts and the theory.

Reprints were made with permission from the publishers.

Contents

1 Introduction

9

2 Non-reflecting boundary conditions (NRBC)

12

2.1 Exact NRBC

12

2.2 Local NRBC

14

3 Absorbing layers

16

3.1 Model problem

18

3.2 Construction of the PML equations

18

3.2.1 First order formulation

19

3.2.2 Perfect matching

20

3.3 Well–posedness of the PML

21

3.4 Stability of the Cauchy PML

22

3.5 Stability of the PML for IBVPs

23

3.6 Stability of the discrete PML

26

4 Finite difference methods for second order systems

29

4.1 The wave equation

30

4.2 Integration–by–parts

31

4.3 Summation–by–parts

31

4.4 Weak boundary treatment with SAT

32

5 Summary of papers

34

5.1 Paper I

34

5.2 Paper II

35

5.3 Paper III

35

5.4 Paper IV

36

5.5 Paper V

36

5.6 Paper VI

37

6 Summary in Swedish

39

References

42

1. Introduction

Wave motion is a dominant feature for problems in many branches of engi- neering and applied sciences. For instance, both industrial materials such as steel pipes and plates [78, 93], and natural structures like the Earth’s surface can support propagating waves that can lead to failure or disaster. Numerical simulations of propagating waves can serve as a complement to theoretical and experimental investigations, and can possibly treat many more important sce- narios that can not be investigated by theory or experiments. Thus, providing valuable information which can be useful in developing modern technologies, exploring natural minerals from the subsurface, and understanding natural dis- asters such as earthquakes and tsunamis.

A defining feature of propagating waves is that they can propagate long dis- tances relative to the characteristic dimension, the wavelength. As an example, strong ground shaking resulting from an earthquake or a nuclear explosion can be recorded far away in another continent sometime after it has occurred. In practice, because of limited resources, computer simulations of such prob- lems are restricted only to areas of interest, where the effects of the strong ground motion is significant. This is typical of numerical simulations of many wave propagation problems including seismic imaging, wireless communica- tion and ground penetrating radar (GPR) technologies.

In numerical simulations, large spatial domains must be truncated to fit into the finite memory of the computer, by introducing artificial boundaries. One can immediately pose the question:

Which boundary conditions ensure that the numerical solution of the initial boundary value problem inside the truncated domain converge to the solution of the original problem in the unbounded domain?

Provided the numerical method used in the interior is consistent and stable,

the numerical solution will converge to the solution of the original problem

in the unbounded domain only if artificial boundaries are closed with ‘accu-

rate’ and reliable boundary conditions. Otherwise, waves traveling out of the

computational domain generate spurious reflections at the boundaries, which

will then travel back into the computational domain and pollute the solution

everywhere. It becomes apparent that the most important feature of artificial

boundary conditions is that, all out–going waves disappear (or are absorbed)

without reflections. Therefore, an effective numerical wave simulator should

be able to treat both physical and artificial boundaries efficiently.

One objective of this thesis is to design effective artificial boundary condi- tions suitable for numerical solutions of partial differential equations (PDE) describing wave phenomena. The effort to design efficient artificial bound- ary conditions began over thirty years ago [35] and has evolved over time to become an entire area of research. Artificial boundary conditions in general can be divided into two main classes: absorbing or non-reflecting boundary conditions (NRBC) [35] and absorbing layers [20]. We also note that for time-harmonic problems there are accurate and efficient artificial boundary procedures, see [107, 47]. Problems of this class are not considered in this thesis.

Most domain truncation schemes for numerical time-dependent wave prop- agation have been developed for constant coefficient wave propagation prob- lems. Many of these methods such as high order local NRBCs and the per- fectly matched layer (PML) which we will discuss in Chapter 2 and Chap- ter 3, respectively, are efficient for certain problems, particularly the scalar wave equation and the Maxwell’s equations in isotropic homogeneous me- dia. For many other problems in this class, such as the linearized magneto–

hydrodynamic (MHD) equations and the equations of linear elasticity, there are yet unresolved problems, see for example [17, 11].

Artificial boundary conditions for more difficult problems such as variable coefficient and non-linear wave propagation problems are less developed. In practice, ad hoc methods are still in use. However, a better understanding of the mathematical properties of the corresponding linear problems will enable the development of efficient boundary conditions for problems of this class.

Second order hyperbolic systems often describe problems where wave prop- agation is dominant. In the numerical treatment of second order systems, the equations are usually rewritten and solved as a first order system, introducing additional degrees of freedom. This is in part due to the maturity of the the- ory and numerical techniques developed in the computational fluid dynamics (CFD) community. Rewriting second order systems to first order system can lead to a loss of computational efficiency [72, 50]. The second objective of this thesis is to derive high order accurate and strictly stable (or time–stable) finite difference approximations for systems of second order hyperbolic partial dif- ferential equations on bounded domains. Numerical methods for second order hyperbolic systems are briefly presented in Chapter 4.

In this thesis, we consider linear time–dependent constant coefficient wave propagation problems in second order formulation. Note that many linear wave equations that appear as first order systems can be rewritten in sec- ond order formulations. Our focus is on the scalar wave equation, curl–curl Maxwell’s equations and the equations of linear elasticity. However, with lim- ited modifications the results obtained in this thesis can be applied to other wave equations.

We begin Chapter 2 by illustrating the fundamental ideas underlying the

derivation of NRBCs and reviewing some well known results of NRBCs. In

Chapter 3, a model problem for second order anisotropic wave equations is

considered, and the PML equations in second order and first order formula-

tions, respectively, are derived. Some of the mathematical properties of the

models are also discussed. Chapter 4 is devoted to a short presentation of

numerical methods for second order hyperbolic systems. A summary of the

included papers is presented in Chapter 5.

2. Non-reflecting boundary conditions (NRBC)

The amount of literature on non-reflecting boundary conditions is enormous and constantly increasing. In this section we attempt a brief review of NRBCs for the wave equation. Elaborate discussions can be found in the review papers [107, 41, 46]. A NRBC, as the name implies, is a boundary condition imposed on an artificial boundary to ensure that no (or little) spurious reflections occur from the boundary [41]. These boundary conditions can be classified into exact (or global) NRBCs, and approximate (or local) NRBCs

. If a boundary condition is such that the artificial boundary appears perfectly transparent, it is called exact. Otherwise it will correspond to a local (NRBC) approximation and generate some spurious reflections.

2.1 Exact NRBC

To begin, we consider the second order scalar wave equation in Cartesian co- ordinates,

1

c

u

= u

+ u

, for t > 0, (x, y) ∈ R

, u = u

, u

= v

, at t = 0.

(2.1)

Here, u is the wave field, c is the wave speed, t denote time, (x, y) are the spatial coordinates, and u

, v

are the initial data. The wave equation (2.1) can describe pressure waves in a homogeneous isotropic media. After some as- sumptions, the wave equation (2.1) can also describe shear waves propagation in elastic solids or wave propagation in an electromagnetic media.

Assume we want to compute the solution of the wave equation (2.1) in the half plane x < 0. In order to close the statement of the problem and make ac- curate computations, a NRBC is needed at the artificial boundary at x = 0.

Exact boundary conditions for (2.1) were first derived in the pioneering work by Engquist and Majda [35] and have recently been reviewed [41, 46] from a modern perspective. The construction of the absorbing boundary conditions

1We note that all exact NRBCs are global, but all global conditions are not exact. Similarly all local NRBCs are approximate but all approximate NRBCs are not local.

[35] uses the fact that any right–going solution u (x, y,t) to (2.1) can be repre- sented by a superposition of plane waves travelling to the right. Such solutions are described by

ˆ

u (x, k

, ω) = a

e

q

(

)

. (2.2)

Here, a

is the amplitude of the wave, and (ω, k

) are the duals of (t, y), sat- isfying ω

/c

− k

> 0 with ω > 0, and the solution ˆ u (x, k

, ω) is the wave function in the Fourier space. By Engquist and Majda [35], the correct bound- ary condition which annihilates all right–going waves at x = 0 is

∂

∂ x u ˆ (x, k

, ω) = −i ω c

s

1 −  k

c ω



ˆ

u (x, k

, ω) , x = 0. (2.3) Notice that the boundary condition (2.3) is prescribed for ˆ u (x, k

, ω) and not for the time–space dependent wave funtion u (x, y,t). But ˆ u (x, k

, ω) is related to u (x, y,t) via the inverse Fourier transform. Therefore, in order to obtain a boundary condition for u (x, y,t) we need to invert the the Fourier transform in (2.3). In theory we can always compute the inverse Fourier transform to determine ∂ u (x, y,t) /∂ x. Unfortunately, the inverse Fourier transform of (2.3) yields the operator

∂

∂ x u (x, y,t) = − F



i ω c

s

1 −  k

c ω



ˆ

u (x, k

, ω)



 , x = 0, with F

u ˆ (x, k

, ω) =

Z

−∞

Z

−∞

e

(

) ˆu(x,k

, ω) dk

dω,

(2.4)

which is nonlocal in both time and space. This is manifested in the difficulty in inverting the pseudo–differential operator √

1 − s

, (with s = k

c/ω), which does not have an explicit local representation. Engquist and Majda acknowl- edged this difficulty and instead resorted to some approximations of √

1 − s

which yield a local differential operator upon inversion. The accuracy of the boundary condition then depends on how well √

1 − s

is approximated. As we will see later, this is the basic idea behind local NRBCs.

The contribution [43, 44, 45] by Grote and Keller is another pioneering work in NRBC. Using the Dirichlet to Neumann (DtN) map, Grote and Keller derived the first exact NRBC on a spherical boundary. We note that the Grote and Keller NRBC is inherently three dimensional, since it is based on special properties of the spherical harmonics. It is believed that no such boundary conditions can be constructed in two space dimensions. This is related to the fact that the Green’s function associated with the wave operator in two space dimensions has an infinitely long tail.

Similarly, starting from the Helmhlotz equation,

ξ

u ˆ = ˆ u

+ ˆ u

, (2.5)

Hagstrom [52] employed the DtN technique to derive the exact boundary con- dition

∂ u

∂ x + 1 c

∂ u

∂ t + F

|k

|

K (|ck

|t) ∗ F u
= 0, x = 0, K(t) ≡ J

(t)

t = 1 π

Z

−1

q

1 − ρ

cos ρtdρ,

(2.6)

which only involves a Fourier transform in the tangential direction and a con- volution in time. It has been reported that by the use of fast algorithms, to- gether with the Fast Fourier Transform, the boundary condition (2.6) can be imposed directly [9].

2.2 Local NRBC

Due to the non-locality of the exact boundary condition (2.3), Engquist and Majda [35] proposed the first hierarchy of local NRBC by replacing √

1 − s

with some rational approximations before inverting the Fourier transform. Ap- proximating √

1 − s

by Taylor expansions and including only the first term yields

1

c u

+ u

= 0, x = 0, t ≥ 0. (2.7) This is the so–called first order Engquist–Majda boundary condition. The boundary condition (2.7) remains exact for a two-dimensional wave propagat- ing normal to the boundary. Including the second term of the expansion yields the second order Engquist–Majda boundary condition

1 c

u

+ 1

c u

− 1

2 u

= 0, x = 0, t ≥ 0. (2.8) We see that including higher order terms of the Taylor expansion to increase the accuracy of the boundary condition in turn introduces higher order deriva- tives at the boundary. However, the inclusion of higher order terms of the Taylor expansion to improve accuracy of the approximation ceases to yield a well-posed problem. This can be cured by the use of Padé approximations, though. The Padé expansion is not the only possible choice. Other expansions like Chebyshev approximations have been studied [106]. Higdon [57] has a more general representation of these boundary conditions

( cos α

c

∂

∂ t + ∂

∂ x ) · · · ( cos α

c

∂

∂ t + ∂

∂ x )u = 0, x = 0, t ≥ 0, (2.9)

where α

· · · α

, are arbitrary parameters. The second order Engquist–Majda

boundary condition (2.8) corresponds to α

= α

= 0 radians. The higher

order derivatives appearing in these boundary conditions greatly complicate

their use in any numerical scheme. As a result, first and second order boundary conditions are most commonly used in practice.

In 1993, Collino [27] introduced a smart idea to remove higher order deriva- tives while retaining high order accuracy at the boundary. The fundamental idea lies in the approximations of √

1 − s

by Padé expansions, and conse- quently introducing a sequence of auxiliary variables φ

. Collino’s idea gave the opportunity to improve the work of Engquist and Majda. Local NBRCs sharing this structure are often referred to as high order local NBRCs, see [41].

Many different high order local NRBCs have been proposed in the past fifteen years, see [41, 53, 48, 107, 46]. The most important property of these high order local NRBCs is that arbitrary order of accuracy can be achieved by in- troducing more auxiliary variables. However, the convergence of the boundary conditions depends on the convergence of the Padé expansion.

Starting from a complete plane wave representation of the time-dependent wave field incorporating both the propagating and evanescent modes, Hagstr- om et al. derived a new local NRBC [55] that annihilates outgoing waves in both propagative mode and evanescent mode regimes. They also presented strong numerical results indicating the long time efficiency of their model.

The general structure shared by all high order local NRBC is that no deriva- tives beyond second order appear, and there are no normal derivatives of any of the auxiliary variables φ

. We point out that on a boundary which has corners, special corner conditions must be used in order for these boundary conditions to yield a well-posed problem.

The major difference between exact NRBCs and approximate NBRCs is

that exact NRBCs are non-local in both time and space. Storage requirement

remains a drawback for exact NBRCs. Another practical challenge in the im-

plementation of exact NBRCs is the shape of the boundary. The boundary

can be very complex such that it becomes extremely difficult to have an ex-

plicit representation of the boundary conditions. We also note that while the

discussion here is formulated in two space dimensions, the ideas carry over

immediately to three space dimensions.

3. Absorbing layers

Another approach to truncate an unbounded spatial domain is to surround the computational (truncated) domain with an artificial absorbing layer of finite thickness, see figure 3.1. All absorbing layers are constructed by modifying the underlying equations such that solutions in the layer decay rapidly. For

Figure 3.1. A computational domain surrounded by the PML.

this method to be effective, it is important that all waves traveling into the layer, independent of angle of incidence and frequency are absorbed with- out reflections. This approach is analogous to the physical treatment of the walls of anechoic chambers. Absorbing layers with these desirable features are called perfectly matched layers (PML). By perfect matching we mean that the interface between the computational domain and the layer exhibits a zero reflection coefficient

. In practice one takes advantage of this property by reducing the width of the layer dramatically, then choosing the damping coef- ficient as strongly as possible to minimize the reflections from the edge of the layer. The accuracy of the entire scheme is then determined by the numerical method used in the interior.

The PML was first introduced for the Maxwell’s equations in the seminal paper [20] by J-P. Bérenger . In its original form [20], the PML is derived by

1A more general interpretation of perfect matching is that the restriction of the solution to the PML problem in the interior coincides with the solution to the original problem, see [10].

splitting the wave functions (magnetic and electric fields) before special lower order terms that simulate the absorption of waves are added. This makes it possible for the PML to absorb all waves exiting the computational domain without reflections, independent of angle of incidence and frequency. There are other formulations of the PML which do not require the unphysical split- ting of the field variables, see for instance [24, 52, 10]. These are the so-called unsplit PML. The unsplit PML is more elegant mathematically. It is also more straightforward to be implemented in a computer program. The PML, both split [20] and unsplit [24, 52, 10], have also been extended to many appli- cation areas like acoustics, elasticity and quantum dynamics. In general, the PML can be interpreted as a complex change of spatial coordinates for the wave equation in the Fourier (or Laplace) space.

The main focus of this thesis is on the PML. However, we also mention that there are other absorbing layers that have been developed, which do not have the perfect matching property. Before the emergence of the PML, Israeli and Orzag [60], Kosloff and Kosloff [66] had begun the construction of ab- sorbing layers for wave equations. In [28, 29], Colonius and Ran proposed an absorbing layer (the super-grid scale model) for compressible flows, by stretching the grid and filtering high frequency components. The paper [12]

by Appelö and Colonius extended the work of Colonius and Ran to linear hy- perbolic systems, in first and second order formulations. In [34], Efraimsson and Kreiss performed a semi-dicrete analysis of a scalar linear model of an absorbing layer (buffer zone). Using their theoretical results Efraimsson and Kreiss suggested how to choose the layer parameters in order to enhance per- formance, then they applied the result to a fully non-linear problem (the Euler equations). Unlike the PML, these absorbing layers do not require the use of auxiliary variables. By construction the super-grid scale models [28, 29, 12]

are linearly stable. They can be used for non-linear problems where the notion of perfect matching is obscure. However, for linear problems, it is doubtful whether these layers can compete with the PML.

In this chapter, we review the PML. To begin, we introduce a simple model

problem for second order strongly hyperbolic systems in two space dimen-

sions. We will then derive a PML model using a plane wave decomposition

that leads to stretching the physical spatial coordinates onto a carefully chosen

complex contour, where spatially oscillating solutions are turned into expo-

nentially decaying solutions. The second order model problem will be re-

written as a first order system. We will also derive a PML model correspond-

ing to this first order system. Finally, we will comment on perfect matching,

well-posedness and stability of the models.

3.1 Model problem

Consider the model problem, the so-called anisotropic scalar wave equation, u

= u

+ u

+ (αu

)

+ (αu

)

, α , x, y ∈ R, |α| < 1, (3.1) describing electromagnetic waves propagating in an anisotropic dielectric me- dia. If |α| < 1, it is easy to show that (3.1) is strongly hyperbolic, see Paper II.

3.2 Construction of the PML equations

The idea is to introduce new coordinates defined by special complex metrics, see [102, 24]. To begin with, we take the Fourier transform in time

ˆ

u (x, y, ω) = Z

−∞

u (x, y,t) e

dt, and consider the time-harmonic problem

− ω

u ˆ = ˆ u

+ ˆ u

+ (α ˆ u

)

+ (α ˆ u

)

. (3.2) The wave equation (3.2) admits plane wave solutions,

ˆ

u (x, y, ω) = u

e

.

Here, (k

, k

) = (k

/ω, k

/ω) are real, u

is a constant amplitude, and (k

, k

) is the wave vector.

Let us consider the case where we want to compute the solution in the half–

plane x ≤ 0, and the PML is introduced outside that half–plane, that is, in x > 0. We now look for the modification of the wave equation such that it has the exponentially decaying solutions

ˆ

v (x, y, ω) = u

e

e

,

in the PML. Here, Γ (x) is a real valued non–negative increasing smooth func- tion, which is zero for x ≤ 0. The decaying solution can be rewritten as

ˆ

v (x, y, ω) = u

e

 k₁

x+^Γ(x)_iω  +k₂y

.

This is a plane wave solution to the wave equation in the transformed coordi- nates ( ˜ x, y), where

˜

x = x + Γ(x) iω .

Next, we analytically continue the wave equation (3.2) onto the complex con- tour ˜ x where spatially oscillating solutions are turned into exponentially de- caying solutions. We have the complex coordinates transformation

∂

∂ ˜ x = 1 s

∂

∂ x , s

:= d ˜ x

dx = 1 + σ

(x)

iω , where σ

(x) = dΓ (x)

dx . (3.3)

In (3.3), σ

(x) ≥ 0 is called the damping function. The PML is then derived by first replacing ˆ u with ˆ v in (3.2), and applying this complex change of variables (3.3) to the wave equation (3.2) (defined for ˆ v), thus yielding

− ω

v ˆ = 1 s

 1 s

v ˆ



x

+ ˆ v

+ 1

s

(α ˆv

)

+

 α 1

s

v ˆ



y

. (3.4)

We note that in order to enhance the absorption and stability properties of the layer more complicated complex metrics s

have been proposed in literature, see [11] and Paper I. Notice that the plane wave satisfying (3.4) is

ˆ

v (x, y, ω) = e

u

. (3.5) Also note that ˆ u solves the wave equation (3.2) in the half–plane x < 0 and ˆ v is the solution of the PML equation (3.4) in x > 0. The wave equation (3.2) is perfectly coupled to (3.4) with the coupling condition

ˆ

u (0, y, s) = ˆ v (0, y, s) , u ˆ

(0, y, s) = ˆ v

(0, y, s) . (3.6) Note that the coupling condition (3.6) is achieved if σ

(0) = 0. We localize the PML in time by first introducing the auxiliary variables,

φ = ˆ 1 s

ˆ v

iω , ψ = ˆ v ˆ

iω .

Inverting the Fourier transforms yields the time–dependent PML equations v

+ σ

v

= v

+ v

+ (αv

)

+ (αv

)

− (σ

φ )

+ (σ

ψ )

, φ

= v

− σ

φ ,

ψ

= v

.

(3.7)

We note in passing that when σ

(x) ≡ 0 we recover the wave equation (3.1).

3.2.1 First order formulation

In order to demonstrate the construction of the PML for a first order system, we introduce extra variables and rewrite (3.1) as a first order system in time and space.

u

= Au

+ Bu

, with u = (u

, u

, u

)

, (3.8) where A =





0 1 0

1 0 0

α 0 0



 , B =





0 0 1

α 0 0

1 0 0



 , and

u

= u

, u

= u

+ αu

, u

= u

+ αu

.

To derive the corresponding PML model for the first order system (3.8) we take the Fourier transform in time and apply the complex change of variables (3.3) in the x-direction. Then choosing the auxiliary variable w and inverting the Fourier transform, we have

v

= Av

+ Bv

− σ

Aw,

w

= v

− σ

w. (3.9)

Also we comment that the auxiliary variables in (3.7) and (3.9) will influence the solutions, v and v, only in the layers where σ

6= 0.

3.2.2 Perfect matching

The most important property of the PMLs (3.7) and (3.9) is the perfect match- ing. This means that the restriction of the solutions to (3.7) and (3.9) in the half–plane x < 0 coincides with the solutions to (3.1) and (3.8), respectively.

There are two standard methods [10, 102] that have been used to study the perfect matching property of the PML. The approach [102] uses plane wave analysis and only accounts for propagating modes. Here, we use the tech- nique [10] which is rooted in the construction of the general solution to the wave equation in the Laplace–Fourier space. The technique [10] is more gen- eral since it includes both the propagating mode regime and the evanescent mode regime. To start with, we take the Laplace tranform in time t → s and the Fourier transform in the tangential direction y → ik

. The equations (3.7) and (3.9) are perfectly matched if

¯

u (x, ik

, s) = ¯ v (x, ik

, s) , and u ¯

(x, ik

, s) = ¯ v

(x, ik

, s) , x ≤ 0,

¯u (x, ik

, s) = ¯v (x, ik

, s) x ≤ 0. (3.10) The variables u, v, u, v are related to ¯ u, ¯ v, ¯u, ¯v via the inverse Laplace–Fourier transformation. Observe that the perfect matching of the second order PML (3.7) requires continuity of the solution and its normal derivative across the interface, while the first order PML (3.9) is perfectly matched if the solutions are continuous across the interface. We can construct modal solutions

¯

u = e

u ¯

(ik

, s) , ¯u = e

¯u

(ik

, s) , (3.11) for the problem in the half–plane x ≤ 0, and

¯

v = e

(

0σ₁(z)dz

) ¯u

(ik

, s) , ¯v = e

(

0σ₁(z)dz

) ¯u

(ik

, s) , (3.12)

for the PML problem in the half–plane x ≥ 0. Direct calculations show that

the first order PML (3.9) is perfectly matched by construction for arbitrary

σ

, while perfect matching is achieved if σ

(0) = 0 in the second order case

(3.7). In computations however, additional smoothness at the interface is often

beneficial.

Remark 1 The modal PMLs [52, 11, 10] are derived by first assuming a modal ansatz of the form (3.12). From the modal ansatz the complex met- ric (3.3) is derived. The PML is then constructed by performing the complex change of variables.

3.3 Well–posedness of the PML

An important mathematical property of a partial differential equation is well–

posedness. By a well–posed problem, we mean that there is a unique solution which depends continuously on the data of the problem. To be precise, con- sider the Cauchy problem

u

= P

 ∂

∂ x , ∂

∂ y



u, t ≥ 0, u(x, y, 0) = u

(x, y).

(3.13)

Here, (x, y) ∈ R

are the spatial coordinates and the symbol P (∂ /∂ x, ∂ /∂ y) denotes the spatial operator. The Cauchy problem (3.13) is weakly (resp.

strongly) well–posed if for every t

≥ 0,

||u (t)||

≤ Ke

||u

||

, (3.14) for u

given in the Sobolev space H

, with s > 0 (resp. s = 0). Here κ and K are independent of u

and t

.

We demonstrate the well–posedness of the second order and first order PML models (3.7) and (3.9). By introducing auxiliary variables we can rewrite (3.7) as a first order system in time and space.

U

= A

U

+ A

U

+ σ

A

U, (3.15)

A

=









0 1 2α 0

1 0 0 0

0 0 0 0

0 0 0 0







 , A

=









0 0 0 1

0 0 0 0

1 0 0 0

1 0 0 0







 , A

=









−1 0 0 0

0 −1 0 0

0 0 0 0

0 0 1 0







 .

It is easy to show that ∀S = (S

, S

) ∈ R

normalized to satisfy

S

+ S

+ 2αS

S

= 1, (3.16) the matrix ˆ A = S

A

+ S

A

has real eigenvalues and a complete system of eigenvectors. It follows that the PML model (3.7) is strongly hyperbolic, thus strongly well–posed, see [49]. However, because of the lower order term σ

A

U the system (3.15) can have solutions that grow in time.

Now consider the PML (3.9) for the first order system (3.8). We can rewrite the PML model (3.9) as follows

V

= B

V

+ B

V

− σ

B

V, (3.17)

B

= A 0 I 0



, B

= B 0 0 0



, B

= A 0 0 0

 .

Since the matrix A in (3.8) has a zero eigenvalue, the matrix ˆ B = S

B

+ S

B

is not diagonalizable. Therefore, the PML model (3.9) is weakly hyperbolic, and weakly well-posed. This is also true for Bérenger’s PML.

If we can rewrite (3.1) as a first order system such that the coefficient matri- ces A and B are invertible it is possible to construct a strongly hyperbolic PML.

This is possible for the acoustic (and advective acoustic) wave equations. But for more complex problems such as the elastic wave equations, rewriting the second order system as a first order system will introduce a non propagating mode which will lead to a weakly well-posed PML, see [17, 11].

However, many examples and computations in literature [19, 26, 11, 17]

indicate that, for linear constant coefficient problems, loss of strong well- posedness may not be ‘disastrous’. We note that if PML models derived for linear constant coefficient problems are to be used for variable coefficient or non-linear problems, it is essential that PML equations are strongly well- posed.

It is also important to note that the discussions here are for Cauchy prob- lems. When physical boundary conditions are important the analysis must be adjusted by considering the theory for IBVPs [74, 75], see also Papers IV and V.

3.4 Stability of the Cauchy PML

Here, we are interested in the temporal behavior of the solution of the PML in the absence of physical boundaries. A typical set–up is a computational do- main completely surrounded by the PML as in figure 3.1. For time–dependent problems, it is not sufficient that the PML is well–posed, it must also be stable.

By definition well–posed problems support exponentially growing solutions.

This is of course undesirable of an absorbing model. The Cauchy problem (3.13) is weakly (resp. strongly) stable if for every t

≥ 0,

||u (t)||

≤ K (1 + γt)

||u

||

, (3.18) for u

giving in the Sobolev space H

, with s > 0 (resp. s = 0). Here K, γ are independent of u

and t

.

A necessary condition for weak stability is that all eigenvalues λ

of the symbol P (ik

, ik

) have non–positive real parts

Perhaps, the earliest (numerical) instability of the PML are reported in [58, 4]. When using the (split-field) PML for the linearized Euler equation, Hu [58] employed a filter to ensure that all fields in the layer decay with time. In [4], Abarbanel and Gottlieb performed a mathematical analysis of Bérenger’s PML for Maxwell’s equation and found that the PML is only weakly well–

posed, and under certain perturbations the solutions of the PML can grow

exponentially in time. Because of the result [4] several unsplit PML were de- veloped [1]. In the paper [3], Abarbanel et al. studied the longtime stability of these unsplit PML. The result is that these layers suffer from long time growth due to the Jordan block present in the lower order damping term. Bécache and Joly [19] used standard Fourier techniques and energy methods to establish the well–posedness and stability of Bérenger’s PML for Maxwell’s equations.

They reported that though Bérenger’s PML is weakly well–posed, the damp- ing term appearing in the layer can at the worst lead to a linear growth. In a subsequent paper [18], Bécache et al. showed that the introduction of the complex frequency shift [77] eliminates the long time growth in the unsplit PML.

While the growth found in [4] is not fatal (since it appear after a very long time), the growth observed in [58] can be destructive. It is inherent in the un- derlying physical problem. Similar growths were also reported for the elastic wave equation in [17]. In [17], Bécache et. al. established an important but negative stability result. They found that the shape of the slowness curve for an arbitrary first order hyperbolic system determines whether a stable split–field PML can be constructed. This is related to the existence of propagating modes with oppositely directed phase and group velocities, that is the so–called back- ward propagating modes. It turns out that this is also true for the modal PML, see [11, 10]. In Paper I, we have also shown that the second order PML for elastic wave equations can support growing solutions when this geometric stability condition is violated.

However, for some (simple) problems such as (3.1), there are linear trans- formations which can modify the slowness diagrams such that a stable PML can be constructed, see [14, 30]. For more complex systems like the elastic wave equations such linear transformations may not be possible. We conclude that the stability properties of the continuous first order PML model (3.9) and the second order PML model (3.7) are linearly equivalent. The PMLs (3.7) and (3.9) can support growing solutions if α 6= 0. As we commented earlier this problem can be cured by introducing extra terms in the PML which modifies the slowness diagrams, see [14] or introducing such a linear transformation first before the PML is derived [30].

3.5 Stability of the PML for IBVPs

For time–dependent PDEs, the introduction of boundary conditions often leads

to a more complicated mathematical problem to be analyzed. In the discrete

setting, more challenging numerical questions can also arise. A PML which

is Cauchy stable can exhibit growth if boundary conditions are included. For

instance in isotropic elastic waveguides, not only that backward propagating

modes can be supported [97]. In the discrete setting, if the boundary condi-

tions are not well treated, the PML can pollute the numerical solutions every- where [97, 12].

For first order hyperbolic PDEs, the theory of IBVPs is rather well–developed.

The Friedrichs theory [37, 38, 39] for symmetric systems with maximally dis- sipative boundary conditions was available, but not sufficient for non-symmetric systems or more general boundary conditions. A more general theory is the Kreiss theory [70, 49, 68] based on the normal mode analysis [7].

The Kreiss technique can be summarized by the following. Consider the problem (3.13) on the upper half–plane (∞ < x < ∞, 0 < y < ∞), with boundary conditions at y = 0,

Lu = g(x,t), y = 0. (3.19)

Here L is a linear boundary operator and g(x,t) is the boundary data.

To begin with,

• localize the problem by freezing all coefficients

• split the problem into a Cauchy problem and a half–plane problem with boundary conditions at y = 0.

The corresponding Cauchy problem can be analyzed using standard Fourier methods. The half–plane problem is analyzed by a Laplace–Fourier technique.

That is, we take the Laplace transform in time t → s, and the Fourier transform in the tangential direction x → ik

, thus yielding an ODE

s ˆu (ik

, y, s) = P

 ik

, d

dy



ˆu (ik

, y, s) + ˆf (ik

, y, s) , ℜs > 0, L ˆu (ik

, 0, s) = ˆg (ik

, s) ,

(3.20)

to be solved on the half–line 0 < y < ∞. The complex function ˆf (ik

, y, s) is given by the initial data. The ODE (3.20) with homogeneous initial ˆf (ik

, y, s) ≡ 0 and boundary ˆg (ik

, s) ≡ 0 data, corresponds to a non–linear eigenvalue problem with s as an eigenvalue. If there are non–trivial solutions of the eigen- value problem with ℜs > 0, the IBVP (3.13) with (3.19) is ill–posed. The boundary condition (3.19) is not useful in a practical sense, and is therefore discarded.

Solving the eigenvalue problem (3.20) with homogeneous data we arrive at the linear system

ϒ (k

, s) δ = 0, ℜs > 0, (3.21) where ϒ (k

, s) is a square matrix determined by the boundary operator L, δ is a set of parameters to be determined. To ensure well–posedness, we require δ ≡ 0 for all ℜs > 0, thus obtaining the determinant condition

F (k

, s) ≡ det ϒ (k

, s) 6= 0, ℜs > 0. (3.22)

The IBVP (3.13) with (3.19) is not well-posed if there are roots s with ℜs > 0 satisfying the characteristic equation

F (k

, s) ≡ det ϒ (k

, s) = 0. (3.23) The Kreiss theory has recently been extended to second order hyperbolic sys- tems [75]. However, the PML is a more complicated system. For second order systems, the PML is a combination of a second order system of PDEs for the transformed wave equation and a first order system of PDEs for the auxiliary differential equations. Unfortunately no stability theory exists for the corre- sponding IBVP including the PML.

Skelton et al. [97], investigated the stability of the PML for elastic waveg- uides, in the frequency domain. They found that an elastic waveguide bounded by free–surface boundary conditions can support backward propagating modes, and these modes deteriorate the performance of the PML. However, they sug- gested that the problem could be cured by constructing a frequency–dependent PML with a damping whose sign depends on the frequency. It is non–trivial to extend their results to the time–domain. In Papers III, IV and V we study the stability of the PML for IBVPs.

In order to demonstrate the stability of the PML for IBVPs, we will again consider the model problem (3.1) on the upper half–plane 0 < y < ∞. At y = 0 we set the flux to zero, having

u

+ αu

= 0, y = 0. (3.24)

This is a model for the free–surface boundary condition in linear elastody- namics. Note that the wave equation (3.1) with the boundary condition (3.24) satisfies a strict energy estimate. However, if we use the physical boundary condition (3.24) with the PML (3.7) the solution explode immediately

, see figures 3.2 and 3.3. In fact, the system (3.7) with (3.24) is unstable in the sense of Kreiss.

Note also that when the physical boundary condition (3.24) only is used in the PML, the derivation of the PML is not consistent on the boundaries. We know that the PML is derived by the complex change of variables ∂ /∂ x → 1/s

∂ /∂ x, defined in equation (3.3). Since the boundary condition (3.24) con- tains u

, in order to ensure a consistent derivation of the PML on the bound- ary we must modify the boundary condition (3.24) by applying the complex change of variables (3.3) to u

. Firstly, we take Fourier transform in time and apply the complex change of variables (3.3) to the physical boundary con- dition (3.24). Secondly, we chose auxiliary variables and invert the Fourier transforms. Depending on the choice of our auxiliary variables we obtain two

2Note that when α 6= 0 the Cauchy PML can support growth independent of the boundary conditions. However, we have used |α| = 0.5 and 10 grid points in the PML such that the discrete Cauchy PML is stable, see Paper II for more details.

possible mathematically equivalent boundary conditions

v

+ αv

+ σ

ψ = 0, y = 0, (3.25) or

v

+ αv

− σ

α φ = 0, y = 0, (3.26) in the PML. When σ

(x) ≡ 0 we recover the physical boundary condition (3.24). Note that both (3.25) and (3.26) are mathematically equivalent. Using the Kreiss technique, we can show that the problem (3.1) with the boundary condition (3.25) or (3.26) does not support temporally growing modes. How- ever, an SBP–SAT, described in the next chapter, approximation of the PML (3.7) with the boundary conditions (3.25), (3.26) results in two different dis- crete systems of equations. Using (3.26) results in temporally growing discrete solution while (3.25) supports no growing solution, see figure 3.2. For more details and results we refer the reader to Papers IV and V.

0 50 100 150 200

10

10

10

10

time

energy

BC1 BC2

Physical BC

Figure 3.2. Time history of the L

norm of the energy kuk with various boundary con- ditions in the PML. The curve BC1 corresponds to using the transforemd boundary condition (3.25), the curve BC2 corresponds to using the transformed boundary con- dition (3.26) and the curve Physical BC corresponds to using the physical boundary condition (3.24) without modification.

3.6 Stability of the discrete PML

In literature, very little attention has been paid to the study of stability of dis-

crete PML models. This is in part due to the fact that if a continuous PML

model supports growth, it may be very difficult to construct accurate and sta-

ble discrete approximations. However, in the discrete setting, derivatives are

time = 1.5

0.05 0.1 0.15 0.2

time = 3

0.05 0.1 0.15 0.2 0.25

time = 1.5

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

time = 3

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

Figure 3.3. Snapshots of the wave field |u| (with |α| = 0.5) in a rectangular waveg- uide bounded on the top and the bottom by the boundary conditions (3.24) in the y–direction. In the x–direction the domain is truncated by a PML containing 10 grid points. The upper panel corresponds to using the transformed boundary condition (3.25) in the PML. The lower panel corresponds to using the physical boundary con- dition (3.24), without the PML transformation. Initially, both solutions are the same.

As time passes the physical boundary condition (3.24) in the PML initiates an explo- sive growth, which propagates first inside the PML, and then spreads everywhere.

approximated, for example by finite differences. In addition, the discrete PML is a finite width layer containing only a few grid points. For a given discretiza- tion the instability or stability in the continuous model can be strengthened or weakened. Therefore a continuous analysis, still very useful, might be incom- plete or too restrictive for the discrete PML.

A relevant article to this study is the recent work [2]. In this paper, Abar- banel et al. performed a systematic experimental study of the long-time behav- ior of discrete unsplit PML models for the Maxwell’s equations. They found that the long-time stability of a given unsplit PML also depends on the chosen discrete approximation.

The study of discrete stability of the PML is also a topic of Papers I, II

and V. In Paper II, we determined the number of grid points needed to ensure

a stable discrete PML for several anisotropic elastic materials violating the

geometric stability condition. In Paper V, we formulate transformed free–

surface boundary conditions in the PML, for elastic wave equations suitable

for the SBP–SAT schemes.

4. Finite difference methods for second order systems

Second order hyperbolic systems often describe problems where wave propa- gation is dominant.Typical examples are the scalar wave equation, Maxwell’s equations, the elastic wave equations and Einstein’s equations of general rel- ativity. In the numerical treatment of second order systems, the equations are often rewritten and are solved as a first order system. This is in part due to the maturity of the theory and numerical techniques developed for CFD.

There are several benefits with solving the equations in second order formu- lation, though. One advantage is that less degrees of freedom are required, and for the same accuracy we need fewer number of grid points per wavelength.

A second advantage is that we can avoid the spurious high frequency modes which can be poisonous when using central difference schemes for first or- der systems. However, while the theory and numerical methods for first order hyperbolic systems are well developed, numerical techniques to solve second order hyperbolic systems are less developed.

For wave propagation problems, it has been demonstrated that high-order accurate time marching methods as well as high-order spatially accurate schemes are more efficient for problems on smooth domains [50, 69]. A major diffi- culty for high order finite difference schemes is the numerical treatment of boundary conditions. For other methods such as finite element methods [61]

and spectral element methods [63, 65], numerical enforcement of boundary conditions can be more straightforward. In addition, the flexibility of numeri- cal computations on unstructured meshes allows for the resolution of compli- cated geometries. However, computational efficiency has continued to make computations on structured meshes attractive. For instance, numerical algo- rithms formulated with finite difference approximations are more intuitive, simple to analyze, and can be readily implemented.

For high order finite difference methods, a methodology that ensures sta-

ble boundary treatment is the SBP-SAT scheme [84, 85], see also Papers III

and VI. The SBP-SAT scheme uses high order finite difference operators that

satisfy the so-called summation-by-parts (SBP) rule. Boundary conditions are

imposed weakly using the Simultaneous Approximation Term (SAT) method

[23]. The scheme leads to a strictly stable approximation. In this chapter, we

will briefly describe the SBP-SAT schemes. For more elaborate discussions

we refer the reader to the suggested references above.

4.1 The wave equation

To begin, we consider the scalar wave equation in one space dimension, u

= u

+ f (x,t), x ∈ (a, b) , t > 0,

u (x, 0) = u

(x) , u

(x, 0) = v

(x) . (4.1) Here, u

(x), v

(x) and f (x,t) are smooth functions which are compactly sup- ported in (a, b) . In order to obtain a well-posed problem, we augment (4.1) with compatible boundary conditions,

α

u

− u

+ γ

u = 0, x = a,

α

u

+ u

+ γ

u = 0, x = b, (4.2) where, α

≥ 0, γ

≥ 0 (i = 1, 2) are real parameters.

To begin with, we introduce the standard L

scalar product and the corre- sponding norm

(u, v) = Z

a

vudx, kuk

= (u, u) . Let us defined the energy by [32]

E

(t) := ku

k

+ ku

k

+ γ

|u (a) |

+ γ

|u (b) |

. (4.3) Multiplying (4.1) with u

, adding the conjugate of the product, and integration–

by–parts yield

d

dt E

(t) ≤ ( f , u

) + (u

, f )

≤ 2ku

kk f (t)k.

We write E

(t) = p

E

(t) p

E

(t). Using the fact p

E

(t) ≥ ku

k, yields p E

(t) ≤ p

E

(0) + Z

0

k f (τ)kdτ. (4.4)

The energy (4.3) is bounded. Thus the problem (4.1) with (4.2) is well–posed.

If we can derive discrete approximations of (4.1) with (4.2), and a correspond- ing discrete energy satisfying (4.4) we say that the discrete approximation is strictly stable. This means that the numerical approximations do not allow any growth not called for by the partial differential equations. For time–dependent problems, strict stability is an important property, particularly if long-time calculations are desired. We note that for IBVPs, it can be difficult to derive strictly stable and higher order accurate schemes.

In order to derive strictly stable discrete approximations, it is important

that the discrete operators approximating the spatial derivatives ∂ /∂ x, ∂

/∂ x

mimic the properties of the continuous operators. Note that the main tool in

deriving the estimate (4.4) is the integration–by–parts property.

4.2 Integration–by–parts

Consider the scalar field u(x,t) defined on the unit interval x ∈ [0, 1]. The spatial operators ∂ /∂ x, ∂

/∂ x

satisfy

 ∂ u

∂ x , u



= Z

0

u ∂ u

∂ x dx {integration–by–parts}

= −

 u, ∂ u

∂ x



+ u(1)

− u(0)

,

(4.5)

 ∂

u

∂ x

, u



= Z

0

u ∂

u

∂ x

dx {integration–by–parts}

= −

∂ u

∂ x

2

+ u(1) ∂ u

∂ x (1) − u(0) ∂ u

∂ x (0).

(4.6)

It is possible to construct discrete approximations of the spatial derivatives

∂ /∂ x, ∂

/∂ x

mimicking the integration–by–parts properties (4.5), (4.6) in a certain discrete norm. Such operators are called summation–by–parts (SBP) operators, and are topics of several works, see for instance [87, 50, 100, 86, 84, 85, 99]. Below, we will briefly describe SBP operators.

4.3 Summation–by–parts

Let u(x,t) be discretized on the uniform grid x

= ( j − 1) h, h = 1

N − 1 , j = 1, · · · , N,

in the unit interval [0, 1]. The grid function at x

is denoted v

(t) and the semi- discrete solution vector is v(t) = [v

(t), v

(t), v

(t), . . . , v

(t)]

. We introduce the discrete scalar product

hv, wi

:= w

Hv, (4.7)

where H is symmetric and positive definite. The scalar product (4.7) defines a norm. Note that with H = hI (where I is the identity operator) we get the standard discrete norm. Let D

, D

be discrete operators approximating the first and the second derivatives respectively, that is

D

≈ ∂

∂ x , D

≈ ∂

∂ x

.

The discrete operators D

, D

are called SBP operators if they satisfy the fol- lowing properties

D

= H

Q, Q

+ Q = B

,

B

= diag(−1, 0, 0, . . . , 1), H = H

, v

Hv > 0 ∀v 6= 0, (4.8)

D

= H

(−M + B

S), M = M

,

v

Mv ≥ 0, H = H

, v

Hv > 0, ∀v 6= 0. (4.9) Here, H defines a norm, Q is almost a skew-symmetric operator, M is called the symmetric part of the operator D

and the operator S is a consistent ap- proximation of the first derivative ∂ /∂ x on the boundaries. Note that D

, D

satisfy

hD

v, vi

= v

HD

v = hQv, vi

= v

−Q

+ B

v

= −hv, D

vi

+ v

− v

, and

hD

v, vi

= v

HD

v = v

(−M + B

S) v

= −v

Mv + v

(Sv)

− v

(Sv)

.

Thus the operators D

, D

completely mimic the integration by part properties (4.5), (4.6) in the H–norm.

Discrete operators satisfying the summation–by–parts properties (4.8) and (4.9) are usually derived using even order standard centered difference schemes in the interior. Close to the boundaries lower order schemes are used such that the operators satisfy (4.8) and (4.9) in a given H–norm. Note that the oper- ator H can be a diagonal norm or a block norm. When using block norms, it can be difficult to prove stability if variable coefficients or curvilinear grids are considered, see [101]. We have used diagonal norms throughout in this thesis.

4.4 Weak boundary treatment with SAT

The construction of SBP operators does not include boundary conditions. For IBVPs, the SBP property alone usually does not ensure a strictly stable scheme.

To ensure strict stability, special boundary treatments are required. Here, we discuss a particular boundary procedure, the SAT method [23]. The main idea is to discretize the PDE (4.1) separately using SBP operators, and the bound- ary conditions (4.2) are discretized using special one–sided difference opera- tors. The semi–discrete PDE and the semi–discrete boundary conditions are then patched together using penalties. The penalty strengths are determined by requiring stability.

The semi–discrete approximation of the boundary conditions (4.2) is B

v

+ B

Sv + B

v = 0, (4.10) where

B

= diag (α

, 0, . . . , α

) , B

= diag (−1, 0, . . . , 1) , B

= diag (γ

, 0, . . . , γ

) ,

and S is the special one–sided approximation of the first derivative ∂ /∂ x, ap- pearing in the second derivative SBP operator D

, defined in (4.9). A SBP–

SAT discretization of (4.1) with the boundary condition (4.2) reads v

=D

v − SAT + F (t) ,

v (0) = u

, v

(0) = v

, (4.11) where the SAT term is defined by

SAT = τ

H

B

v

+ B

Sv + B

v
, and τ

= 1 is a penalty. We can rewrite (4.11) as

v

= Bv

+ Dv + F(t),

v(0) = u

, v

(0) = v

. (4.12) Here,

D = −H

M + B

, B = −H

B

. Let the semi-discrete energy be defined by,

E

(t) := kv

k

+ v

Mv + γ

v

+ γ

v

. (4.13) Multiplying (4.12) with v

H, adding the transpose of the product leads to

d

dt E

(t) ≤ hf, v

i

+ hv

, fi

≤ 2kv

k

kfk

. As before, we write E

(t) = p

E

(t) p

E

(t), and use p

E

(t) ≥ kv

k

, to ob- tain

p E

(t) ≤ p

E

(0) + Z

0

kf (τ) k

dτ. (4.14) The semi–discrete approximation (4.11) is strictly stable.

In Paper VI, we derived SBP–SAT schemes for more general second or-

der hyperbolic systems. For Dirichlet boundary conditions, it is possible to

impose the boundary conditions strongly, by injection. The injection method

for Dirichlet conditions in the SBP framework, for systems of second order

hyperbolic PDEs is also discussed in Paper VI.

5. Summary of papers

This chapter presents the summary of the individual papers included in this thesis. The material contained in this thesis discuss the construction, analysis and numerical implementation of the perfectly matched layer for second order wave equations. We investigate the stability of both the continuous and the discrete PML models. We also derive strictly stable finite difference approx- imations suitable for systems of second order hyperbolic PDEs. The sections below briefly describe each paper, and highlight the main results in the papers.

5.1 Paper I

In this paper, we construct a general PML for second order hyperbolic sys- tems, focusing on the linear anisotropic elastic wave equations in two space dimensions. The layer equations are derived by applying a complex coordi- nate stretching directly to the second order equations. The resulting system is strongly hyperbolic. When the geometric stability condition is violated, by a standard perturbation argument our PML at constant coefficients suffers from the same high frequency instability affecting the modal PML and the split field PML for first order systems. However, in computations using standard second order finite differences, our PML behaves much better than a standard first order PML. We have found several reasons for this.

In a discrete setting the unstable modes may be of higher frequency than can be represented, or well represented, on the grid. The temporal behavior of such modes cannot be expected to be predicted by continuous analysis.

We also show that this effect can be enhanced by coordinate compression and the complex frequency shift in the layer. However, coordinate compression increases the stiffness of the problem.

Secondly, we observe that the geometric instability gives rise to growing modes, with bulk localized to part of the layer, and propagating tangentially.

If a Cartesian domain is surrounded by layers, the bulk of the unstable mode

eventually moves into a corner region, and decays. We analyze the stability

properties of the corner region, as before. The result is that, the corner region

does not support temporally growing solutions. All waves propagating into

the corner region will decay immediately.