ORNUTM-13373

?
On the Parallelization of the
Acoustic Wave Equation with
Absorbing Boundary
Conditions

C. T. White
V. A. Protopopescu
J. Barhen

1
This report has been reproduced directly from the best available copy.

Available to DOE iind DOE contractors from the Office of Scientific and
Technical Information, P.O. Box 62, Oak Ridge, TN 37831; prices available
from (615) 576-8401.

Available to the public from the National Technical Information Service, U S .

Department of Comrnerce, 5285 Port Royal Rd., Springfield, VA 22161.

This report was prepared as an account of work sponsored by an agency of

the United States Government. Neither the United States nor any agency
thereof, nor any of their employees, makes any warranty, express or implied,
or assumes any legal liability or responsibility for the accuracy, completeness,
or usefulness of any information, apparatus, product, or process disclosed, or
represents that its use would not infringe privately owned rights. Reference
herein to any specific commercial product, process, or service by trade name,
trademark, manufacturer, or otherwise, does not necessarily constitute or
imply its endorsement, recommendation, or favoring by the United States
Government or any agency thereof. The views and opinions of authors
expressed herein do not necessarily state or reflect those of the United States
Government or any agency thereof.
 DISCLAIMER

Portions of this document may be illegible

in electronic image products. Images are
produced from the best available original
document.
 ORNL/TM-13373

ON THE PARALLELIZATION OF THE ACOUSTIC WAVE EQU

WITH ABSORBING BOUNDARY CONDITIONS

C. T. White
Department of Mathematics
California Institute of Technology
1200 East California Boulevard
Pasadena, CA 9 1125
E-mail: [email protected]

V. A. Protopopescu
J. Barhen
Center for Engineering Systems Advanced Research (CESAR)
Oak Ridge National Laboratory, P. 0. Box 2008
Oak Ridge, TN 37831-6355
E-mail: vvp @ornl.gov,barhenj @ornl.gov

July 1998

Research supported by the Office of

Fossil Energy, U. S. Department of
Energy under the Advanced Computing
Technology Initiative.
I

Prepared by the
Oak Ridge National Laboratory
Oak Ridge, Tennessee 3783 1
managed by
Lockheed Martin Energy Research Corp.
for the
U. S. DEPARTMENT OF ENERGY
Under Contract No. DE-ACO5-96OR22464
Contents

ABSTRACT 1

1 INTRODUCTION 1

2 BACKGROUND 2
2.1 ABSORBING BOUNDARY CONDITIONS FOR THE
ACOUSTIC WAVE EQUATION . . . . . . . . . . . . . . ............... 2
2.2 THESPONGEENHANCEMENT . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.3 THE IFP-3D CODE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

3 DISCRETE ABSORBING BOUNDARY CONDITIONS 7

4 PARALLELIZATION OF ACOUSTIC WAVE EQUATION 9

4.1 EXACT DIAGONALIZATION OF THE ACOUSTIC WAVE EQUATION WITH
DINCHLET BOUNDARY CONDITIONS . . . . . . . . . . . . . . . . . . . . . . . 9
4.2 IMPOSSIBILITY OF EXACT DIAGONALIZATION OF THE ACOUSTIC WAVE
EQUATION WITH ABSORBING BOUNDARY CONDITIONS . . . . . . . . . . . 11

5 DISCUSSION 14

ACKNOWLEDGEMENTS 17

REFERENCES 19

iii
 ABSTRACT
Many practical problems involve wave propagation through atmosphere, oceans, or terrestrial crust.
Modeling and analysis of these problems is usually done in (semi)infinite domains, but numerical
calculations obviously impose restriction to finite domains. To mimic the actual behavior in the
(semi)infinite medium, artificial absorbing boundary conditions are imposed at the boundaries,
whereby waves can only exit, but not enter the finite computational domain. Efficient absorbing
boundary conditions are difficult to analyze and costly to run. In particular, it is of interest
to assess whether the wave equation with (approximate or exact) absorbing boundary conditions
admits a suitable diagonalization. This would open the possibility for parallelizing many important
numerical codes used in applications. In this paper we propose a set of stable, local, absorbing
boundary conditions for the discrete acoustic wave equation. We show that the acoustic wave
equation with absorbing boundary conditions cannot be exactly diagonalized.

1 INTRODUCTION

Propagation of acoustic, elastic, or electromagnetic waves in oceans, atmosphere, or terrestrial crust

is often cast in the form of semi-infinite or infinite problems. Typical instances occur in underwater
acoustics, seismic wave propagation, and electromagnetic wave generation by cosmic objects or
antennas. However, for computational purposes these practically infinite problems must be recast
within finite domains.
Several approaches have been investigated to address this requirement. One possibility is to map
the (semi) infinite domain into a finite one. Another approach is to use smaller computational
domains by substituting the true equations with simpler ones in the exterior domain. Then one
can use the form of the solution of the simplified equations to construct boundary conditi~nsfor
the true equations in the inner domain.
The most widely used approach has been to introduce artificial computational boundaries while
maintaining the governing differential equations unchanged. On the artificial boundaries one should
design boundary conditions that ensure that the incident waves are fully transmitted and the re-
flection is, ideally, zero. These artificial boundary conditions imposed on the solution of the wave
equation in a (semi)infinite domain are called absorbing boundary conditions (ABCs). Various ABCs
have been used in the numerical modeling of wave propagation (see for instance, [l-151 and refer-
ences therein). Theoretically, pseudodifferential operators can be constructed that ensure perfect
transmission (Le., no reflection) of the outgoing waves at the artificial boundaries. Unfortunately,
except in the one-dimensional case, these boundary conditions are non-local in both space and
time and therefore not practical for numerical calculations. Resorting to local approximations is
somewhat of a gamble since it may lead to ill-posed problems for the wave equation [1,3,7]. Even
when the local approximations lead to well-posed problems, the perfectly absorbing property of the
boundary is limited to certain angles of incidence (usually normal incidence) [1,7,9-121. To enhance
9 the (imperfect) absorptive quality of the artificial boundaries, additional damping and filtering
mechanisms have been devised [4,5,14-161.However, although they reduce reflection - sometimes
to very low levels - dampers and sponges (or filters) introduce an extra computational burden
and may inherently prevent the implementation of more efficient methods, such as parallelization.
Even in the absence of sponges one faces an additional problem. Assuming that local stable ap-
proximations of ABCs have been constructed, these continuous boundary conditions have to be

1
again approximated by discrete ABCs which are then coupled to a standard interior discretization
of the wave equation. Since in t,his process some of the required properties may be lost, it seems
more natural to look directly for local discrete versions of ABCs that lead to stable, efficient, and
possibly parallelizable schemes.
The aim of this paper is to analyze the parallelizability of the acoustic wave equation (AWE) with
ABCs. The paper is structured as follows: In Section 2 we review the AWE and various ABCs
and absorption enhancing techniques. We present a new local, stable, discrete ABC for the exact
discrete AWE in Section 3 and compare its performance with the three dimensional code developed
by the Insitut Franqais du Petrclle (IFP-3D) [5]. Section 4 is devoted to an analysis of the ABCs
as to their potential for parallelisation. The results are analyzed in the final section.

2 BACKGROUND

2.1 ABSORBING BOUNDARY CONDITIONS FOR THE

ACOUSTIC WAVE EQUATION

Let u(z,t)solve the AWE

utt - c2v- vu = g ( 2 , t )
in the domain s1. Usually this equation is supplemented with initial conditions, u(z,t = 0) = uo(z),
ut(z,t = 0) = u1(z) and with boundary conditions at the boundary of R, dR. The typical boundary
conditions for the AWE are Dirichlet, Neumann, and Robin. The homogeneous Dirichlet condition
is imposed at a fixed boundary and is expressed by

ulan = O .
Waves striking a fixed boundary will reflect with an amplitude of opposite sign as the incoming
wave. The homogeneous NeumaIin condition is imposed at a free boundary and is expressed by

where n is the unit vector n0rm.d to the boundary surface. Waves reflect from such a boundary
with the same amplitude as the incoming waves. The Robin condition is given by

where again, n is the unit vector normal to the boundary surface. This condition is also known as
the mixed or elastic boundary condition and in order to obtain stable solutions, X is required to be
greater than zero.
These boundary conditions, are c,bviously non-absorbing. To effectively mimic an infinite physical
domain within a finite computational domain we need artifical boundaries which almost completely
“transmit” waves out of the region of computation. An exact solution to this problem exists only for

2
 the one-dimensional (1D) AWE. In higher dimensions, various approaches have been proposed to
deal with this problem, often based on implementing approximations to pseudodifferential operators
on the boundaries. We give a brief summary of the techniques used in Ref. 5:
The problem domain, R, is assumed parallelepipedic. As such, conditions for the faces, edges, and
corners are required. The homogeneous three dimensional (3D)AWE reads

*
[-- a2 - ( " + - + - ) a2] u = o .aa29
1
c2 at2 0x2 ay2 (5)

To derive a boundary condition for the faces x = X m i n or x = Xm,, we solve formally for 2:

For the face x = Xmin, we select the positive root; for the face x = Xm,, we choose the negative
root. We show the details for the face x = Xmin.
The pseudodifferential operator in the RHS of Eq. (6) is approximated by using the Taylor expansion
for the square root, d G = 1- ;a + O(a2), yielding

For the face x = Xmin, we obtain

typos with signs, they should

be opposite, since for Xmin,
we select positive root
This equation is called second order (15") purmiul approximation of the acoustic waves equation in
the direction x > 0 [5].
For the surface x = Xm,, the corresponding ABC reads:

here is also typo with signs,

they should be opposite.
Further, in (8a), the signs are correct
The ABCs at the other surfaces have similar expressions, obtained by circular permutation of the
coordinates.
Let us now consider the line defined by the intersection of the two surfaces x = Xmin and y = Ymin.
The associated well-posed condition along this line is obtained by satisfying the ABCs used along
each surface and the AWE.
The condition for the surfaces x = Xmin and y = Y- are

3
respectively.
Solving Eq. (8a) for @ and Eq. (8b) for 3,
and substituting into the AWE, Eq. ( 5 ) , we obtain

Similarly, we deduce the conditicns for the others edges.

Finally, let us consider the four corners at the bottom of the model. For example, the ABC at the
corner defined by the intersectior. of the surfaces x = Xmin, y = Ymin, and z = Zmin, is obtained by
the combination of the three AB13 along these surfaces and the AWE. The condition for the lines
{X = Xmin , y = Y,n}, {X = Xmin 7 = & i n } , { y = Y,n = &in} are

respectively. Solving Eq. (loa) fcirs,

Eq. (lob) for @, and Eq. (1Oc) for 3,
and substituting
these expressions into the AWE, Eq. ( 5 ) , we obtain the ABC at the corner {x = Xmin , y =
Ymin , z = Zmin}

No additional conditions are needed for the four corners of the free surface.

2.2 THE SPONGE ENHANCEMENT

An alternative approach to reducing reflection from the boundaries is based on sponge filters that
are located near the houndary of the domain. Their intended effect is to damp waves traveling
from the boundary into the interior (thereby reducing the amplitude of reflected waves), and leave
unaffected the waves approaching the boundary. Of course, the sponge filter may be used in
conjunction with some form of ABCs for enhanced overall effect. For illustration, we derive the
sponge filter for the 1D situation, The 1D AWE can be written as

4
Each of these first order operators in the L.H.S. of Eq. (12)
.~governs the propagation of waves in
2
.

a single direction: waves traveling right satisfy f - = 0, while waves traveling left satisfy
-;g+g=0.
-
Near the right boundary we introduce an artificial source that propagates waves to the left with
a negative amplitude, thus partially or - in the ideal situation - completely canceling the waves
reflected from the boundary. We express this as follows:

----- ; zau
1 a2u a2u - f - € ( x ) ( 1 +z) .
c2 at2 ax2

The term ~ ( xallows

) us to control the damping factor at each spatial point, so that we may take it
to be zero deep in the interior of the domain and increasing as we approach the right boundary. We
note that large changes in E over a short distance tend to introduce additional reflections. Thus,
for the sponge term to be effective, the damping coefficient must increase gradually, often over a
large number of grid points.
In 3D the situation is more complicated since the factorization above does not apply. The goal,
however, remains to damp waves traveling in from the boundaries. For plane geometry, large
surfaces and homogeneous media, the 1D approximation of the sponge at the 5 = X ,, boundary,

---
1 a2u (-+-+s)
a2u = f - € ( x ) (1eaua + $ ) .
c2 a t 2 ax2 ay2

works reasonably well. Similar terms are used for the other faces.
An enhancement to the “one-way” sponge filter described above is obtained by damping both
incoming and outgoing waves near the boundary [4]. While providing a greater net damping effect,
(as waves are damped both before and after they reflect from a boundary), this approach may
require one to disregard data obtained from points within the sponge filter, since even outgoing
waves have been disturbed at these points. .
We construct the “two-way” sponge by simply introducing, as sources, canceling factors for both
incoming and outgoing waves. For the 3D AWE near the x = X,, boundary, we find

1 a2u
- --
c2 a t 2
(-d2U
ax2
+-
ay2
+ -)
a2u a2u
a22
=f - E(Z) (-1 -
c at +
au au
-€(X) (;1 au - s>
au

=f - ;2€ ( X ) - au .
at
The damping term involves now a coefficient that depends on the speed and on the distance to
the boundary. Similar terms describe “two-way” sponges for the other faces.

2.3 THE IFP-3D CODE

We include here some brief comments on the three dimensional code developed by the Insitut
F’ranqais du Petrole (IFP-3D) to solve the AWE for use in geophysical exploration and analysis

5
[5]. This code is designed for seismic studies of subsurface exploration. This feature alone requires
massive calculations to account jor each point, in the discrete 3D domain. Additional features of
the code also play a part in its total computational complexity which results in a rather costly
program.
The main features of the IFP-3D code are:

0 full 3D model of the AWE in a parallelepipedic domain;

0 implements the second order paraxial boundary conditions on five of the six faces of the
domain, implements either the Dirichlet or Neumann condition on the sixth face to simulate
either a fixed or free surfact:;
0 includes an additional one-way sponge enhancement to further reduce reflections;
0 assumes a uniform spatial grid, i.e. A, = Av = A,;
e uses an explicit scheme: a finite difference scheme approximates time derivatives to second
order and space derivatives to tenth order (where possible);

0 is sequential in time: solves updates for all spatial points at a given time step before proceeding
to the next;
e allows the use of several so’urceterms, including a point source given by a Ricker function;
and
0 writes output for specific lines or planes of points in the domain which are the “simulated
receivers.”

The most consequential drawbacks of the IFP code are:

e the explicit method is only conditionally stable and requires rather small At (hence a large
number of time steps) while allowing only relatively low source frequencies;
0 the algorithm is essentially sequential in time: current attempts at parallelizing the IFP-3D
code involve stepping sequentially through time, with costly node communications required
at each step; and
0 there is a high computatioiial cost associated with the code’s approximation of ABCs: in
addition to the often mentilmed cost involved in the paraxial boundary conditions, there is
also the high cost of performing the sponge update on several grid layers (often 20 or more)
near the boundary.

Our research has been motivated in large part by the need to enhance the efficiency of the 3D model
algorithms (in particular the IFP-3D code) while not sacrificing their accuracy. This requires the
development of (hopefully unconditionally) stable algorithms that have potential to support parallel
implementations. The following analysis focuses on these issues.

6
 3 DISCRETE ABSORBING BOUNDARY CONDITIONS

In this section we derive discrete ABCs which are local in both time and space and considerably
simpler in expression than the paraxial equations used in Ref. 5. As mentioned before, our motiva-
tion in constructing these boundary conditions comes from not only the computational complexity
of the paraxial equations, but also from the difficulty encountered when trying to incorporate the
paraxial equations into a time parallel algorithm. Moreover, we have chosen to start directly with
discrete ABCs since, eventually, the numerical solution of the problem requires discretization. We
first give a brief account of the derivation of these conditions, followed by an analysis of their
effectiveness.
In constructing the discrete ABCs, we assume that the 3D (parallelipipedic) domain is discretized
in the form:

(0,1,. .. ,N X + 1) x (0,1,. .. ,N Y + 1)x (0,1,. ..,N Z + 1)

The boundaries, in the form of surfaces, edges, or corners have at least one coordinate equal to
+ +
either 0, N X 1, N Y 1, or N Z 1. +
The new ABCs generalize to the multidimensional case the method of characteristics that works
exactly in the 1D problem. We will now show details of the new boundary conditions for the 3D
problem; a similar procedure applies in the 2D case.
We denote by ubnI(i,j,k) the solution of the discrete 3D AWE at discrete time n and grid point
(i,j,k). The discrete homogeneous 3D AWE obtained from the finite difference method with second
order approximations reads:

where p ( i , j , k) = +,
C2(ij,k)A2
c(i,j , k) is the speed of sound at the grid point (i,j,k),At is the discrete
time step, and h = Az = Ag = A, is the discrete spatial step.
We write ufU+l]and dn-l] by using the ansatz

(i,j , k)
u[n+l] . p ( i , j , k ) [ ~ [ ~ I ( i + l , j , k 1) +(u[nl(i,j+l,k)+u[n](z,j-
= 2 1,k)

+uin1(i,j,k + 1)+ u["](i,j,k - l))]+ (1- 3 p ( i , j ,k)) uIn1(i7j,k )

uIn-1l(i,j,k) + 21 (u["l(i,j+ 1,k) +u["l(i,j - 1,k)

= p ( i , j , k ) [ d q i - 1,j)k)

+dnI(i,j,k+ l)+u[nl(z,j,k- l))]+ (1-3p(i,j,k))uW(i,j,k) .

7
which clearly satisfy Eq. (16).
For i = 0, 1 5 j 5 N Y , 1 5 k 5 N Z , Eq. (17a) becomes

1
+ 1 , k ) +u["](O,j - 1 , k )
u[n+'l(O,j,k) = p(O,j,k)[u[nI(l,j,k) ~ ( u [ " I ( O , j +
+ dn1(D,j,k + 1) + d n ] ( O , j ,k - l ) ) +] ( 1 - 3p(O,j,k))u[n](O,j,k). t

For each face we perform a similar splitting, which selects only the interior point in the direction
normal to the surface. To derilre the edge condtions, we split the two spatial second derivates
normal to the intersecting faces, solving the equation which contains only those interior points. For
example, at i = 0, j = 0, 1 5 k 2; N Z , we obtain

l , k ) + uin](O,
u [ ~ + ' ] ( O , 0, k ) = p(O,O, k ) l k [ n ] ( 0, 1,k) (19)
+ 1
+ +
,(u["](O,0, k 1 ) uin](O,0, k - l ) ) ]+ ( 1 - 3p(O, 0, k ) ) ~ [ ~ l0,
( Ok ), .

Finally, for corner points we split all three of the spatial derivates, again solving the subequation
which contains only these interior points. For the corner i = 0, j = 0, k = 0, we get

+ u [ (0,1,0) +
u[~+'l(o,o,o)= p(o,o,o)[u[~~(l,o,o) (0, 0, l)]
+
( 1 - 3p(0,0,0))u[n](0,0,0).

These new ABCs are.now local in both time and space and are compatible with the discretized
equation used for the interior of the domain.
To analyze the effectiveness of these ABCs we compared them with the second order paraxial
conditions used in the IFP-3D algorithm and with the Dirichlet condition. We used two discrete
domains with 150 x 150 x 100 and 150 x 150 x 80 grid points, respectively. The source is located
(75,75,75) in the center of the domain in the X Y plane, 25 (or 5) points from the bottom face
boundary. A single receiver is placed initially at the same point as the source; on later shots it is
moved along the 2 = 75 plane t o realize various angles of incident waves namely 0", 17", 25", and
41". The speed of sound, c, is equal to 750.00m/s, At = .167 x lO-*s, A, = Ay = A, = 3.0rn, and
the central frequency of the Ricker source is 20.00Hz.
In evaluating the effectiveness OF the new ABCs, we considered the stated objectives: simplicity,
speed, and potential for paralleljzation. Specifically, one must decide if the absorption of the new
ABCs reaches an acceptable level, and if the greater absorption provided by the paraxial conditions
is not offset by their complexity In particular, for the serial implementation we are interested in
the time necessary to calculate the boundary update at each time step. Also, since the complex
paraxial ABCs cannot be easily incorporated into a fuzzy parallel scheme, one must assess whether
the simpler ABCs will allow thc simultaneous diagonalization of the resulting system, providing
full space and time parallelism.
Figures 1 and 2 show the relative reduction in amplitudes of reflected waves for both the paraxial
and new ABCs. One sees a greater reduction from the paraxial conditions, but both significantly

8
reduce the unabsorbed Dirichlet reflection. Since the tests were performed on relatively small
scale models, we have no definitive conclusions regarding the computation times of the various
boundary conditions. Preliminary data show that the new ABCs are calculated slightly faster than
the paraxial ABCs. Due to locality, this effect would be significantly enhanced by parallelization.

4 PARALLELIZATION OF ACOUSTIC WAVE EQUATION

The implementation of time parallelism for numerical solutions of partial differential equations has
been shown to hold great potential [19-211. Whenever applicable, this formalism would enable the
use of fully implicit methods which, due to their unconditional stability, dramatically reduce the
number of required time steps and provide a massive degree of coarse-grain temporal parallelism
with minimal communication and synchronization requirements. Such algorithms have proven to
be highly suitable for implementation on massively parallel MIMD architectures such as the Cray
T3D and the Intel Paragon [19-21 and references therein]. We shall assess the applicability of this
methodology to the AWE with ABCs, with specific implementation targeted at the enhancement of
the IFP-3D code. The approach we follow is to search for the eigenvalue-eigenvector decomposition
of the matrices resulting from the finite difference method. To illustrate this procedure we use the
AWE in a finite uniform domain with Dirichlet boundary conditions and then turn to the situation
involving ABCs.

4.1 EXACT DIAGONALIZATION OF THE ACOUSTIC WAVE EQUATION

WITH DIFUCHLET BOUNDARY CONDITIONS

Once again, assume that the parallelipedic domain is discretized in the form

(0,1,. . .,NX + 1) x (0,1,...,N Y + 1) x (0,1,.. ., N Z + 1)

Recall that the boundaries are comprised of those points having at least one coordinate equal to 0,
+ +
N X 1, N Y 1, or N Z 1. +
We consider the discrete system obtained through the finite difference method,

where

1. is an ( N X ) x ( N Y ) x ( N Z ) - vector (note that we may omit boundary points from the

vector u since we are assuming Dirichlet (zero) boundary conditions);
2. L is the block tridiagonal [ ( N X )x ( N Y ) x ( N 2 ) ] 2discrete Laplacian matrix [19];

+
3. p = czAz is constant throughout the domain, and h = A, = Ay = A,; and

4. fin] is the source function at time t = nAC\t.

9
 Time steps

. Fig. 1. Amplitude of reflected waves as a function of time for the 150 x 150 x 80 domain at 0":
Dirichlet (l),IFP3D (2), new ABCs (3).

Time steps

Fig. 2. Amplitude ofreflected waves as a function of time for the 150 x 150 x 100 domain at 0":
Dirichlet (1), IFP-3D (Z), new ABCs (3).

10
The Crank-Nicholson approximation:

results in the system

The eigenvalue-eigenvector decomposition of the discrete Laplacian L is known analytically [19].

We do mention that L is diagonalized by an orthogonal transformation matrix 8 resulting from a
discrete sine transform. Let 8LQ-l = A be diagonal, where of course 8-1= Q T .
Therefore,

By making the substitution ] , obtain the decoupled system

= @ u [ ~we

where = [email protected] that the eigenvalues of L are negative, so that ( I - &A) is nonsingular.
Let D = ( I - $PA)-' and rewrite the above system as

If the above tw&erm recursion is homogenous, its iteration involves a linear combination (with
coefficients which depend upon the initial conditions) of powers of the matrices D zk (D2-
which may be computed very efficiently since D is diagonal. If, on the other hand, the above is
a nonhomogeneous recursion, we may still efficiently apply standard algorithms such as recuvsive
doubling or cyclic reduction, since D is diagonal. All that remains is to use the stored matrix
8-1 = QT to perform the transformation of v back to u. Finally, note the unconditional stability
of this method resulting from the Crank-Nicholson approach.

4.2 IMPOSSIBILITY OF EXACT DIAGONALIZATION OF THE ACOUS-

TIC WAVE EQUATION WITH ABSORBING BOUNDARY CONDITIONS

The diagonalization approach presented in the previous section can be carried out for periodic
and Neumann boundary condition, but cannot be immediately generalized to more complicated
boundary conditions. Indeed, the difficulty in diagonalizingthe AWE with the second order paraxial
boundary conditions was among the motivating factors of our present analysis. We shall show
that, the paraxial boundary conditions as well as our other proposed schemes fail to meet specific
conditions necessary in order to diagonalize exactly the corresponding AWE problem either in
contunious or discrete form. We now turn to an analysis of these conditions.

11
First we consider an AWE which is continuous in time and discrete in space. I n general, this can
be written in the form:

a2u du
X- -t. Y- = 2u, for X,Y ,2 E MN(R) ,
at2 at
where MN(R) denotes the set of the N x N matrices over the real numbers $2.
We require that Eq. (21) with suitable ABCs satisfy the following three conditions:

I. closely approximates an undamped discrete AWE in the interior of the domain; this ensures
that the ABCs do not distort the interior solution.

11. provides absorbing, or nori-reflecting boundaries, i.e. approximates well a (semi)infinite

medium; and
111. admits a full, simultaneous diagonalization of the three matrices involved which ensures the
potential for parallelization. The simultaneous diagonalization requires the existence of two
matrices 8, 9 E MN(R),9 invertible, such that OXQ,8YS,82Q are all diagonal matrices.

Before establishing conditions necessary for a triple Y ,2) to meet these three criteria, we state
(X,
a weaker version of I1 which will be useful later.

11’. for any initial condition, Eq.(l) and g(z, t ) = 0 has a unique solution, u(t),such that limt-oo 11
u(t) I]= 0. This conditions ensures that after a finite time, the wave will have essentially left
the finite medium (no reflections).

Suppose that 11’ and I11 hold for some X,Y ,2 E MN(R). Let 8,Q be the matrices from I11 above
and define QXQ = a, 8YKP = I?, @2Q = A.
Defining v(t) = W1u(t),we obtztin the following decoupled system:

dv
+ -a2v
-
at2
+ r - =at~ v .
The differential equation associated to a typical mode is of the form
q5V” + yv‘ = xu .

Since Eq. (21) is second order in t [me,Eq. (23) will also be specified with two initial data. Therefore,
in order to ensure the existence arid uniqueness of solutions to Eq. (23)) and hence Eq. (21), it follows
that # # 0.
Indeed, if # = y = X = 0, then amy twice differentiable function v ( t )satisfymg the initial conditions
also satisfies the (trivial) DE. Tlius, solutions are not unique. On the other hand, if q5 = 0 while
not both y and A are zero then che DE in Eq. (23) is of order one or zero. Since two initial data
have been specified, this problem. is in general overdetermined and solutions need not exist.
From this observation, we conclude that both X and 8 are nonsingular. Without loss of generality,
therefore, we may assume X = I .

12
 We now consider solutions to Eq. (23) where 4 # 0. Since

we obtain the series of implications

for each mode (24)

This last implication enables us to conclude y # 0 and X # 0. Indeed, if y = X = 0, then solutions

to Eq. (23) are linear in time and in general do not decay. If y # 0, X = 0, then solutions to
Eq. (23) involve a constant plus an exponential; certain choices of initial conditions guarantee that
this constant is non-zero. Finally, if y = 0, X # 0, the roots of the characteristic equation of
Eq. (23) are A m . If Re (m) # 0, then one of the two exponentials will grow unboundedly
(m)
with time. On the other hand, if Re = 0, then the solution will have a constant magnitude,
which will be nonzero for some initial conditions.
Thus, in every case except y # 0 and X # 0, the mode in question fails to decay in time for certain
initial conditions. Then, by the above remark, 11 u(t)11 does not converge to zero. Thus, we must
have y # 0 and X # 0; hence r and A, (consequently 7 and 2 ) ,are non-singular.
We now argue that 7 and 2 commute, and thus we may take Q = @-I. Indeed, assuming the
simultaneous diagonalization as before, we obtain

Finally, since any finite collection of pairwise commuting, diagonalizable matrices can be simulta-
neously diagonalized via a similarity transformation, we may take Q = 8-l.
In summary:
If we suppose that the system

d2U du
-+y-=2u
at2 at
is (2) is simultaneously diagonalizable and (ii) has solutions that vanish for large t , then the following
necessarily hold:

1. y and 2 are non-singular;

2. y2 = 2 y ; and
a
3. The system can be simultaneously diagonalized by a similarity transformation.

As before, 2 is the discrete Laplacian matrix for the interior points (plus some suitable expression
for the boundary points), and hence it cannot be a diagonal matrix. Therefore, since y commutes

13
with 2, if y is a diagonal matrix then y = T I N for some y f irz. But this will introduce a constant
absorption of the waves throughcut the whole medium and not only near the boundaries. Thus for
y # 0, condition I will not be satisfied.
We turn now to an AWE that is discretized not only in space but also in time and is represented
by the second order recursion relation

for A , B , C E MN(R).
As before, we suppose there exist 8,Q f M N ( % ) ,@ invertible, such that d = @A@, = 8B@,
C = OCQ are all diagonal matrices. The resulting decoupled system is

where = @I-'u['I. Each modc is then a second order recursion of the form

where 210 and v1 are specified initial conditions.

The analysis proceeds as in the continuous case. If a = p = y = 0, then solutions to this trivial
recursion relation are not unique. If either Q = 0, p # 0, y # 0 or y = 0, Q # 0, /3 # 0, then the
above is actually a first order recurrence. The problem is overspecified since two initial terms are
given, and in general, solutions need not exist. The same applies if two of the three coefficients are
zero. Thus, we conclude that ths matrices d and C are nonsingular. Therefore A, C, and 8 are
also nonsingular.
Without loss of generality we can take A = I . Furthermore, in order that 8 and @ exist, it is
necessary that B and C commute, hence I , 23, and C may be simultaneously diagonalized via a
similarity transformation. Therefore, we may assume Q = 8-l.
The general solutions of the recursion relation Eq. (27) with Q = 1 is

where 7-1, r2 are the roots of the characteristic polynomial x2 + pa: + y = 0, and A1 and A2 are
complex constants determined bjr the initial conditions.
Thus, in order that -+ 0 as n -+ 00, irrespective of the initial conditions, a final necessary
condition is

In summary, three conditions similar to (1) - (3) in the continuous case are shown to be necessary
here as well. Also, as in the previous case, this will contradict at least one of the requirements I -
111.

14
 5 DISCUSSION

We conclude with a short discussion of these necessary conditions and, in particular, their impact on
exact diagonalization of the AWE with ABCs. We shall first demonstrate how the ABCs discussed
in this paper fail to meet at least one of the necessary conditions, thus establishing the impossibility
of simultaneous diagonalization and parallelization for these ABCs.
The second order paraxial boundary conditions fail to meet the nonsingularity condition. Note, for
example, that the corner condition involves only the first and second time derivatives of u but not
the function itself; therefore, in the notation of the preceding section, the matrix 2 contains some
rows comprised entirely of zeros and hence is singular.
The general sponge filter fails to meet the commutativity condition. Since a sponge filter involves
only multiples of a time derivative at each point, the matrix y (in the notation of the previous
section) is diagonal. As it is that the matrix 2 performs the discrete Laplacian updates on interior
points, 2 cannot itself be diagonal. Thus, y and 2 commute if and only if y = T I N for some
constant 7,and of course we are only interested in y > 0. Thus, the only sponges (used indepen-
dently of other boundary routines) which satisfy the necessary conditions are those which have a
constant damping factor. (Note that the constant factor throughout the domain is the coefficient
c: kA2
of the time derivative, -+&jk , and not simply E . )
Finally, the new ABCs fail to meet the nonsingularity condition as discussed in the previous section
for the time-discretized system. Specifically, this condition requires that the matrices operating on
the vectors uin+l]and u [ ~ -be~ ]nonsingular. Since the new ABCs make use of only the two time
+
steps (n l)A, and nAt in computing each update, the matrix operating on the vector uln-'1
contains many rows consisting entirely of zeros and hence is singular.
As a means by which to introduce the general problem, consider first the special case of a uniform
sponge, or, using an earlier notation, a constant diagonal matrix Y . This method is an unacceptable
solution to our problem, as it fails to restrict the region of damping to points close to the boundary,
but instead damps waves uniformly throughout the domain. Thus, any damping effect one wishes
to create on or near the boundary will also be created throughout the entire domain. Consequently,
if the damping coefficient is large enough to reduce reflections to an acceptable level, it will similarly
damp interior propagation and thereby destroy the accuracy of the interior solution.
The analysis of a general sponge (and matrix y ) is an open problem that should be addressed after
criteria I - I1 are given a more precise quantitative statement.
We also leave open at this point, the complete conclusions of our analysis of second order diagonal-
izable systems for solving the AWE. Whether or not there are any workable schemes which satisfy
our necessary conditions is not known nor guessed at this point. One item to notice is that in the
case of a diagonal y , the requirements that y be nonsingular and commute with 2 are somewhat
competing notions, forcing the condition of uniformity which results in either insufficient damping
at the boundary or excessive damping in the interior. We suggest that a similar incompatibility may
exist for the more general situation. However, it is our belief that an analysis of this last question
will rest upon a clear and quantitative reformulation of the requirement that any scheme respect
the accuracy of the original AWE in the interior of the domain. We recommend this direction for
future research.

15
E
 ACKNOWLEDGEMENTS

This research was performed at the Center for Engineering Science Advanced research, Oak
Ridge National Laboratory. Funding was provided by the Office of Fossil Energy under the
Advanced Computing Technology Initiative under contract DE-AC05-960FC22464 with
Lockheed Martin Energy Research Corporation. The authors thank Drs. D. B. Semeraro and D.
B. Reister for their help, and Dr. Leigh House (LANL) for continued support.

17
 REFERENCES
1. B. Engquist and A. Majda, Absorbing Boundary Conditions for the Numerical Simulations
b
of Waves, Math. Comp, 31, 629 (1977).
2. B. Engquist and A. Majda, Radiation Boundary Conditions for Acoustic and Elastic Wave
Calculations, Comm. Pure and Appl. Math, 32, 313 (1979).
3. M. Israeli and S.A. Orszag, Approximation of Radiation Boundary Conditions, J. Comp.
Phys., 41, 115 (1981).
4. C. Cerjan, D. Kosloff, R. Kosloff, and M. Reshef, A Nonreflecting Boundary Condition for
Discrete Boundary Condition for Discrete Acoustic and Elastic Wave Equations, Geophysics,
50, 705 (1985).
5. R.L. Higdon, Absorbing Boundary Conditions for Difference Approximations to the Multi-
dimensional Wave Equation, Math. Comp., 47, 437 (1986).
6. R.L. Higdon, Numerical Absorbing Boundary Conditions for the Wave Equation, Math.
Comp. 49, 65 (1987).
7. J. Sochacki, R. Kubichek, J. George, W.R. Fletcher, and S. Smithson, Absorbing Boundary
Conditions and Surface Waves, Geophysics, 52, 60 (1987).
8. A. Bamberger, P. Joly, and Y.E. Roberts, Second Order Absorbing Boundary Conditions for
the Wave Equations: A Solution for the Corner Problem, SIAM J. Numer. Anal., 27, 323
(1990).
9. L.T. Long and J.S. Liow, A Transparent Boundary for Finite Difference Wave Equation,
Geophysics, 55, 201 (1990).
10. R.L. Higdon, Radiation Boundary Conditions for Elastic Wave Propagation, SIAM J. Numer.
Anal., 27, 831 (1990).
11. D.R. Burns and R.A. Stephen, Three-Dimensional Numerical Modeling of Geoacoustical Scat-
tering from Seafloor Topography, J. Acoust. SOC.Am., 88, 2338 (1990).
12. D. Givoli, Non-Reflecting Boundary Conditions, J. Comp. Phys. 94, 1 (1991).
13. L. Anne, J. Brac, and H. "ran, Propagation 3D d'Ondes Acoustiques, Institut Fkanqais du
Pdtrole, April 1993.
14. Chengbin Peng and M. Nafi Toksoz, An Optimal Absorbing Boundary Condition for Finite
Difference Modeling of Acoustic and Elastic Wave Propagation, J. Acoust. SOC. Am., 95,
733 (1994).
15. J.P. Berenger, A Perfectly Matched Layer for the Absorption of Electromagnetic Waves, J.
1
Comp. Phys., 114, 185 (1994).
16. J. Nordstrom, Accurate Solutions of the Navier-Stokes Equations Despite Unknown Outflow
Boundary Data, J. Comp. Phys., 120, 184 (1995).
17. P. Joly and J. Tuomela, A New Theoretical Approach to Absorbing Layers, SIAM J. Numer.
Anal., 34, 671 (1997).

19
18. Xiaobing Feng, Absorbing Boundary Conditions for Electromagnetic Wave Propagation,
Math. Comp., to appear.
19. N. Toomarian, A. Fijany, and J. Barhen, Time-Parallel Solution of Linear Partial Differ-
ential Equations on the Intel Touchstone Delta Supercomputer, Concurrency: Practice and 4

Experience, 6 , 641, (1994).

20. A. Fijany and P. C. Messina, Unconditionally Stable Explicit Method for Massively Parallel I
Solution of Acoustic Wave 'Equations, SPIE, 2571, 212 (1995)

21. A. Fijany, M. A. Jensen, Y. M m a t - S a d , and J. Barhen, A Massively Parallel Computation

Strategy for FDTD: Time and Space Parallelism Applied to Electromagnetics Problems, IEEE
Trans. Ant. Prop. 43, 1441 (1995).

I
 INTERNAL DISTRIBUTION

1-5. J. Barhen 25. N. S. Rao

6. Y . Braiman 26. D. B. Reister
7. G. A. Geist 27. M. S . Smith
8. C. W. Glover 28. C. Touzet
9. W. C. Grimmell 29. M. A. Unseren
10. H. K. Liu 30. S. Wiseman
11. E. M. Oblow 31. T. Zacharia
12. L. E. Parker 32. Central Research Library
13. N. Peterfreund 33. ORNL Laboratory Records - RC
14 - 24. V. A. Protopopescu 34-35. ORNL Laboratory Records - OSTI

EXTERNAL DISTRIBUTION

36. Dr. Robert E. Price, ER-15, Office of Basic Energy Sciences, Department of Energy,
Washington, DC 20585
37. Dr. Fred Aminzadeh, UNOCAL, 14141 SW Freeway, Suite 301-225, Sugarland TX 77478.
38. Dr. John Blair, JBX Technologies, 25 Moore Road, Wayland, MA 01778
39. Amir Fijany, Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove
Drive., Pasadena, CA 91 109
40. Dr. Oscar P. Manley, 13212 Wye Oak Drive, Gaithersburg, MD 20878.
41. Dr. Charles R. Weisbin, Robotics and Mars Exploration Technology Program, Jet Propulsion
Laboratory, 4800 Oak Grove Drive; Pasadena, CA 9 1 109
42 - 47. C. White, Mathematics Department, California Institute of Technology, 1200 East California
Blvd. Pasadena, CA 91 125