\documentstyle[11pt,times]{cmu-art}
\input{psadobe}
%\input psfig
\input epsf
\pagestyle{plain}
\include{ogmacro}
\pagestyle{plain}
\setcounter{topnumber}{3}
\renewcommand{\topfraction}{1}
\setcounter{bottomnumber}{3}
\renewcommand{\bottomfraction}{1}
\setcounter{totalnumber}{3}
\renewcommand{\textfraction}{0}

%\setlength{\textheight}{9.0in}
%\setlength{\textwidth}{6.5in}
%\setlength{\footheight}{0.0in}
%\setlength{\topmargin}{0.0in}
%\setlength{\headheight}{0.0in}
%\setlength{\headsep}{0.0in}
%\setlength{\oddsidemargin}{0in}
%\setlength{\parindent}{1.5pc}
%\setlength{\flushbotton}

\begin{document}
\bibliographystyle{plain}
\include{psfig}

%mention something about utah people \cite{sikorski90}
%mention //ellpack!

\begin{center}
{\bf EARTHQUAKE GROUND MOTION MODELING IN LARGE BASINS}
\end{center}

The main objective of the proposed research is to develop the
capability for predicting, by computer simulation, the ground motion
of large basins during strong earthquakes, and to use this capability
to study the seismic response of the Greater Los Angeles Basin. This
problem represents a grand challenge both because of its intrinsic
importance to hazard mitigation, and due to its modeling and
computational complexity.  Anticipated increases in raw computing
power over the next five years will not be sufficient for successful
understanding and prediction of earthquake ground motion, unless
simultaneous advances are made in physical modeling, numerical
approximations, parallel algorithms, and software tools for the new
generation of high-performance parallel systems. Our effort will focus
on three components necessary to achieve this goal:

\begin{itemize}

\item {\em Physical and mathematical modeling:} A basic obstacle is
characterization of the geomorphology of large basins. Considerable
work has been conducted to date, and the main geometric features and
mechanical properties of a number of large basins are known
sufficiently. Our challenge is to incorporate this knowledge into
three-dimensional, heterogeneous models, a task which has only recently
been possible for portions of the largest basins.  Furthermore,
existing models are incapable of predicting any inelastic behavior.
Although computational complexity will increase, the development of
nonlinear models may be necessary in sedimentary regions that
experience strong shaking.

\item {\em Numerical approximation:} Existing models
cannot predict the response of large sedimentary basins because their
use of structured grids necessitates much finer spatial and temporal
resolution than can be accommodated. Unstructured mesh methods for
solving the partial differential equations (PDEs) of wave propagation
must be developed and implemented.


\item {\em Computation:} The full power of the new generation of 
parallel systems has not been harnessed by geophysicists and
earthquake engineers for real problems. To exploit their power, new
parallel algorithms and software tools must be developed. Problem
solving environments must be created to allow scientists and engineers
to solve PDE problems on parallel computers without requiring them to
know the underlying architecture and communication system. Parallel
systems and problem solving environments must be able to handle the
irregular communication patterns that arise from unstructured mesh
methods.
\end{itemize}


A concerted, unified effort is required along all three lines to be
able to predict earthquake ground motion for the largest of seismic
regions, a task that requires teraflop speed and terabyte primary
memory. To this end, a collaboration was initiated in late 1991
between engineers and computer scientists at Carnegie Mellon
University to address unstructured mesh problems and applications. The
proposed research seeks to build upon this initial work in
mathematical models, numerical methods, and software tools, and upon
geophysical research over the last two decades directed specifically
at elucidating the effects of site conditions, path propagation, and
source effects on seismic ground motion. A team comprised of
engineers, geophysicists, computer scientists, and a mathematician,
will attack the problem of predicting strong ground motion in highly
seismic, large regions. During the development and initial calibration
stages, we will concentrate on the Valley of Chusal, a
well-characterized and instrumented small valley located in an
actively seismic region in the Confederation of Independent States (CIS).
Our major goal is to conduct a
numerical simulation of the earthquake response of the Greater Los
Angeles Basin, with a view towards developing more rational
recommendations for prescribing design ground motion, and in
particular, for taking spatial variations into consideration in
seismic design provisions.

It is important to point out that while the proposed work is motivated
by our interest in gaining a better understanding of strong seismic
motion in large basins, from a mathematical point a view the
simulations will entail the numerical approximation of PDEs using
unstructured meshes on parallel computers. The development of parallel
unstructured PDE solvers is of central interest in scientific
computing. Thus, the numerical algorithms and software tools we
develop will be applicable to a wider range of problems that require
unstructured meshes, such as transonic flow, reservoir simulation,
materials processing, and pollutant transport.

\section{Physical and mathematical modeling issues}

Central to the mission of the National Earthquake Hazard Reduction
Program is the reduction of the earthquake hazard to the general
population. To this end, it is essential that within seismic regions
new facilities be designed to resist earthquakes and existing
structures be retrofitted as necessary.  Assessing the free-field
ground motion to which a structure will be exposed during its lifetime
is a critical first step in the design process. Ground motion is
usually specified through seismic design spectra, which essentially
prescribe an equivalent lateral force that the structure must
withstand without failure. This force, which is a function of the
fundamental natural period of the structure, is based upon past
seismic history in the general geographic location and position with
respect to possible earthquake sources such as active faults; expected
earthquake magnitudes; and general geologic conditions. Observations
of ground motion during recent strong earthquakes have shown, however,
that local site effects, which are normally given only passing
attention in design, can be extremely significant, and can adversely
affect structural safety. Three common effects often observed in
basins or sedimentary valleys are an amplification and significantly
longer duration of the surface ground motion with respect to that on
rock, as well as a more rapid spatial variation of the ground motion
which can cause large differential base motion of extended structures
such as bridges or dams.

Examples of all these effects are plentiful. Perhaps the most dramatic
recent occurrences are those in Mexico City in 1985
\cite{Seed88,Sanchez-Sesma88} and within the San Francisco Bay area 
during the 1989 Loma Prieta earthquake \cite{Boore89}. For both of
these events and their aftershocks amplifications greater than 4 or 5
and durations of up to 15 to 30 seconds greater than the corresponding
motion on rock were quite common, due to local site conditions.  From a
detailed study of aftershocks of the 1989 Loma Prieta earthquake in
the Santa Clara valley, Frankel et al.\ \cite{Frankel91} concluded that
the surface waves generated at the valley confluence were largely
responsible for the long duration of shaking observed in the
aftershock seismograms. A similar conclusion involving also the
subbasement was reached by Bard et al.\ \cite{Bard88} and Kawase and
Aki \cite{Kawase89} regarding the 1985 Mexico earthquake. Analogous
studies have been conducted for portions of the Kanto Basin for a
number of several well-instrumented earthquakes, with similar
observations \cite{motosaka92,phillips93}. A recent study was able to
deal with the three-dimensional seismic response of a portion of the
Greater Los Angeles Basin by using a simplified geometry and and a
high velocity acoustic approximation \cite{Graves92}. The study
indicates that the presence of the large sediment-filled basins
amplifies significantly the strength of the waves observed in these
areas, with a complex pattern of shadow zones and bright spots created
in areas both inside and outside of the basins.

A simple illustration of these effects is shown in Figure \ref{xli}b
which depicts the horizontal component of surface displacement in the
longitudinal and transverse directions in an alluvial, homogeneous,
elastic valley of nearly semi-ellipsoidal shape, embedded in a stiffer
valley surrounded by a halfspace, as illustrated in Figure \ref{xli}a.
The excitation is given by an incident plane SV-wave propagating at an
angle of $90^\circ$ with respect to the $x$--axis (North) and an angle
of $20^\circ$ degrees with respect to the vertical axis. The four
synthetic seismograms show, from top to bottom, the prescribed
free-field displacement, the displacement in the E--W direction at
point A within the small valley, the one-dimensional (1D) displacement
for the corresponding flat-layered system, and the transverse
displacement in the N--S direction at the same location. These
results, which we obtained on the CM-2 during the
course of an exploratory study sponsored by NSF \cite{li92},
correspond to a somewhat shallow valley (semiaxes ratios 2:1:0.2) for
which 1D theory might be expected to provide satisfactory results near
the center of the valley\footnote{The normalized physical parameters
in this example are: $\frac{\beta_{1}}{\beta_r}=0.1$,
$\frac{\beta_2}{\beta_r}=0.4$; $\nu_1 = \nu_2 = 0.25$, $\nu_r = 0.3$;
and $\frac{\rho_1}{\rho_r}=0.69$, $\frac{\rho_2}{\rho_r}=0.88$.
$\beta$,$\nu$, and $\rho$ denote shear wave velocity, Poisson's ratio,
and density, respectively, and the subscripts 1 and 2 refer to the
soft and stiffer valleys, and $r$ to the surrounding rock. These
values are intended to represent an inner valley with a medium
impedance contrast with respect to the surrounding larger valley, and
a high contrast with respect to the rock.}.

By comparing the three top displacements it is clear that the presence
of the layers by themselves causes an increase in the amplitude and
duration of the response with respect to the free-field motion. Most
important, however, are the differences between the response of the
flat-layered system and that near the valley center; the strong
additional increase in the response amplitude and duration is due, in
part, to the reverberation of the body waves through the valley, but
mostly due to surface Rayleigh and Love waves generated at the
confluence of the soft valley with the surrounding stiffer media. The
bottom trace in Figure \ref{xli}b arises entirely due to the
three-dimensional character of the response; since the excitation has
components only in the E--W and vertical directions, the response in
the N--S direction, in contrast to the significant motion that takes
place within the valley, would vanish identically for a flat-layered
system.  Figure \ref{xli}c shows the distribution of the displacement in
the E--W direction within the entire valley at a given instant, while
Figure~\ref{xli}d depicts the peak value of this displacement over the
entire duration of the seismic motion within the basin. The former
figure is a snapshot of an animation created to help visualize the
spatial distribution of seismic ground motion over time (see
accompanying videotape). These two figures illustrate dramatically
that large differential motions can take place over a short distance.
As already indicated, this has practical implications for the
design of long structures.

%(The physical parameters in this example, $\frac{\beta_s}{\beta_r} =
%0.33$, $\frac{\rho_s}{\rho_r} = 0.93$, $\nu_s = 0.25$, and $\nu_r =
%0.3$, where $\beta$, $\rho$, and $\nu$ denote shear wave velocity,
%density, and Poisson's ratio, and the subscripts $s$ and $r$ refer to
%the valley and the surrounding rock, are intended to represent a
%valley with a medium impedance contrast with respect to the
%surrounding halfspace).

Most of the computational seismological studies carried out to date on
the response of sedimentary deposits have been made under the
assumption that the soil remains elastic throughout an earthquake.
Some observations of ground response lend validity to this assumption
\cite{rogers84}. There is clear evidence, however, that during
the 1989 Loma Prieta earthquake inelastic soil behavior occurred in
many locations. This behavior was observed primarily but not
exclusively in areas covered by Bay mud or artificial fill, and
involved both liquefaction
\cite{Chameau91} as well as soil amplification \cite{Darragh91,Chin91}.
There is evidence of nonlinear soil amplification
also from earthquakes in Japan \cite{kamiyama92}, involving both soft
and stiffer soil deposits.  An interesting situation is provided by
Chusal Valley, a small, well-instrumented sedimentary valley
consisting of competent soil and located in a seismically active area
in the CIS. The shear-wave velocity within the valley increases
monotonically with depth from 250 m/sec at the surface to 1000 m/sec
at the basement. A number of earthquakes with maximum surface
acceleration ranging from $10^{-5}$g to $10^{-3}$g and from 0.01g to
0.2g have been recorded within the valley. Figure \ref{chusal} shows the
Fourier amplitude ratio of the response in the E--W direction at two
sites within the valley, one in the middle and the other at the edge,
as a function of frequency.  While on the whole the two sets of
spectral ratios corresponding to the weak and stronger earthquakes are
similar, it is noteworthy that the average amplitude of the first
peaks is smaller for the stronger earthquakes than for the weaker
ones. Most noticeable is the reduction in the value of the frequency
at which these peaks occur for the case of the stronger earthquakes.
Since both of these features are characteristic of nonlinear response,
this suggests that Chusal Valley might have responded inelastically to
the stronger set of earthquakes.  While seismologists have only rarely
(e.g. \cite{Joyner75,Joyner75a,yu92}) considered inelastic behavior,
within the geotechnical engineering community it has been generally
accepted for a number of years that soils behave inelastically if the
excitation is sufficiently strong to cause shear strains exceeding,
say, 0.1 percent. This typically can occur for soft to medium stiff
clays and sands if the maximum acceleration on rock exceeds 0.2g
\cite{Seed69}.  Analyses of site response conducted by geotechnical
engineers, however, have been limited to one-dimensional flat-layered
models subjected to plane vertically incident SH-waves. See Finn
\cite{Finn91} and Yu et al.\ \cite{yu92} for state-of-the-art reviews from
a geotechnical engineer's and from seismologists' point of view.

As exemplified by the cases cited in the previous paragraphs it is now
generally recognized that while one- and two-dimensional local models
can help explain observed behavior in certain situations, a complete
quantitative understanding of strong ground motion requires a
simultaneous consideration of effects of earthquake source,
propagation path, and local site conditions
\cite{Campillo89,Chin91,Graves92}. In the case of soft to medium
soils, or if the excitation is unusually strong, consideration of the
inelastic behavior of the soil material also becomes necessary.

Recent advances in the development of constitutive laws for soil
materials make it now possible to take into consideration the
interaction of the pore fluid (water) with the solid soil skeleton.
This two-phase behavior allows the realistic modeling of the evolution
of soil deformation up to the point of failure under transient
excitation. For instance, Zienkiewicz and his collaborators
\cite{Zien90a,Zien90b} have succeeded in performing an
approximate reconstruction of the failure of the lower San Fernando
dam during the 1971 earthquake. They have also validated results of
their computations by comparison with experimental results obtained on
physical models tested in a centrifuge. Such models, however, have
never been used for predicting ground motion in seismic areas
susceptible to site amplification or liquefaction.

\section{Numerical Approximation Issues}

The propagation of seismic waves through the earth is most commonly
modeled by the elastic wave equation, a system of hyperbolic partial
differential equations. Because analytical solution is rarely
possible, the solution must be approximated numerically. See
\cite{aki88,sanchezsesma87,vidale88} for 
surveys of common time domain and frequency domain approaches with
applications primarily to two-dimensional situations. An elastic
finite-difference modeling of the 1971 San Fernando, California
Earthquake, done directly in the time domain, is reported in [..]. The
finite element method approach is illustrated, e.g., in [..].  More
recently, the finite difference method has been used \cite{frankel92}
to simulate three-dimensional seismic waves directly in the time
domain in the Santa Clara Valley from a Loma Prieta aftershock.

Boundary integral equation methods
\cite{cole78,sanchezsesma93,hisada93}, the Aki-Larner wavenumber
method and its modifications \cite{ohori92,uebayashi92}, and other
frequency domain methods (e.g.\ \cite{toshinawa92}) have been used
successfully for purely linear problems of moderate size. The main
advantage of these methods is that they require that the solution be
expressed only in terms of unknown variables defined on boundary
surfaces or interfaces. The difficulty with these approaches is that
they all lead to full systems of equations, which require at least
$O(n^{2})$ operations per frequency for $n$ variables. In addition,
the evaluation of the necessary Green's functions can be tedious.
Thus, their use has thus far been limited to two-dimensional and
low-frequency three-dimensional cases, generally involving simple
structure.  Various methods have been proposed recently in an effort
to make the computations more efficient. On one hand, Gauss
beam-boundary element methods \cite{cerveny82,benites89,kato93} and
full-space Green's functions have been used to simplify the evaluation
of Green's functions [..], while on the other hand localization
techniques have been developed to produce sparse matrices
\cite{rokhlin90,zeng92}. The latter techniques have not yet been used
in seismic motion studies.  Combined domain finite element and
boundary integral methods have also been used
\cite{bielak91,mossessian92} to model basins with inhomogeneous
properties.

To design a numerical approximation scheme capable of dealing with
large, highly heterogeneous basins, several features must be taken
into consideration. These characteristics are shared by many problems
in computational science:

\noindent {\bf Ability to resolve complex geometry.} The domains of many
realistic physical problems modeled by partial differential equations
do not admit simple description. The interface between a sedimentary
valley and the surrounding rock is typically complex, as can be seen
from Figure \ref{kanto}. The ideal numerical scheme should be
able to reproduce this geometry with minimum error. 

\noindent {\bf Ability to resolve multiscale phenomena.} Many physical
problems display behavior on widely differing temporal and spatial
scales, usually due to variations in geometry, coefficients, forcing
functions, or boundary conditions. The difference in material
properties that a wave experiences as it propagates through a soft
sedimentary valley results in spatial resolution requirements for
sediment and surrounding rock that can differ by a factor of 20. The
approximation scheme must be able to adapt to this changing scale of
the physics. We illustrate this idea with two meshes of a
cross-section through a portion of the LA Basin (Figure
\ref{lasection}).  The mesh of Figure~\ref{struc} is typical of finite
difference grids, and is termed ``structured'', in the sense that it
is regular, and thus can be described by a small number of parameters.
The two dashed lines indicate the interface between the two layers of
the valley and the surrounding rock. Although an attempt is made to
tailor the grid spacing to the physics of the problem by stretching
it, it still results in many more equations than desired in the rock
regions.  As a result of its flexibility, the finite element mesh can
be tailored to the shear wave velocity in each layer, as illustrated
in Figure~\ref{unstruc}.  This results in a considerable reduction in
the number of variables.

\noindent {\bf Ability to resolve nonlinear constitutive laws.}
Recent recognition of the possible importance of nonlinearities in
earthquake-induced ground motion (summarized above) calls for the
creation of methodologies to simulate inelastic behavior resulting
from earthquake-induced ground motion. A greater challenge is the
simulation of the highly nonlinear phenomenon of liquefaction.

\noindent {\bf Adaptivity.} Our discussion has emphasized static
meshes---those which remain fixed throughout temporal integration.
Such an approach is sufficient for many problems, since high
resolution areas are known ahead of time, e.g. at reentrant corners
and interfaces. For other problems, however, either such areas are not
known a priori, as in the formation of shocks or slip surfaces, or
these high gradient regions propagate in time and must be tracked. In
such cases, greater serial efficiency can be obtained by adaptively
refining the mesh based on some local or global error estimate. The
best way to perform this is unclear for distributed memory parallel
machines. Enriching the approximation in a subdomain may lead to load
imbalance; attempting to balance the load by redistributing nodes can
be communication intensive. In any case, for linear and mildly
nonlinear waves, a static mesh will suffice; for liquefaction models
and other highly nonlinear physical phenomena, some form of adaptivity
may be required.

In view of these desired features, it appears that unstructured finite
elements can best provide the necessary adaptivity and multiscale
resolution needed to solve truly large problems, a conclusion that has
long been accepted in the field of solid mechanics, and the importance
of which is becoming recognized in many diverse areas\footnote{See,
for example, the National Research Council report on research
directions in computational mechanics \cite{nrc91}}.  Nonlinearity and
inhomogeneity considerations rule out integral equations and other
dimensional reduction methods; irregularity, multiscale behavior, and
complex geometry rule out structured finite differences, spectral
methods, and ``fast'' solvers. Considerations of simplicity and
modularity lead us to unstructured finite element methods for the more
general problems.

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\section{Computational issues}


Realistic ground motion modeling problems often exceed the capacities
of existing supercomputers.  Consider for example the target of our
work, the Greater Los Angeles Basin, depicted in Figure~\ref{LA}. For
a shear wave velocity of 0.4km/s and a period of 0.5s, a finite
difference approximation of the elasticity operator would place grid
points 0.02km apart to achieve second order accuracy. The region of
interest has dimensions 140km$\times$100km$\times$20km; thus, a
regular discretization, governed by the softest layer, requires 35
billion grid points with three displacement components per node. At
least a terabyte of primary memory would be needed, and on the order
of $10^{13}$ operations would be required at each time step. The
stability condition associated with standard explicit time integration
of the semidiscrete equations of motion imposes a time increment at
least as small as 0.004s. Thus, a hypothetical teraflop computer would
have to run for three days to simulate a minute of shaking.  Of
course, the development of unstructured discretization methods and
efficient parallel implicit solvers will result in considerable
reduction. By tailoring mesh size to wavelength, the unstructured mesh
methods we propose will result in a considerable reduction in number
of unknowns, and implicit solvers will allow longer time steps. Even
then, however, the full basin requires the capabilities of
next-generation supercomputers.

In order to model the ground motion in large basins, we must address
two main issues: (1) the lack of parallel algorithms and tools for
unstructured PDE problems, and (2) the need for better communication
systems in the parallel computers that will ultimately run the
earthquake codes.  Addressing these issues will require close
cooperation between users, compiler writers, and parallel system
developers.  Here is a brief survey of the problems we face:

\noindent {\bf Mesh generation.} Unstructured mesh generators for
arbitrary regions are difficult to design.  A good survey of mesh
triangulation methods is provided by Bern and Eppstein~\cite{Bern}.
It is not difficult to automatically triangulate an arbitrary
two-dimensional region.  It is more difficult to find a
tetrahedralization for an arbitrary three-dimensional region.  The
real difficulty, though, is that a mesh generator should avoid
producing elements with large aspect ratios, because angles close to
$180^\circ$ lead to large discretization error~\cite{Babuska}, and
very small angles cause the stiffness matrix to be
ill-conditioned~\cite{Carey}.  A good mesh generator will also
minimize the number of nodes and elements in the mesh.

\noindent {\bf Mesh partitioning, placement and routing.} We partition the
domain into subdomains by breaking the mesh into pieces, one per
processor.  Partitioning unstructured meshes is deceptively difficult:
the pieces should have nearly the same number of elements, ensuring a
good load balance; and each piece should have a relatively small
boundary, minimizing the parallel communication.  Most effective
partitioning algorithms are heuristic in nature (e.g.
\cite{PSL90,byun92,farhat88}), and have superlinear behavior---they
become painfully slow as mesh size increases. Partitioning must be 
repeated often for moving meshes, making a fast partitioner
especially important for such problems.

\noindent After partitioning, each subdomain is mapped to a separate
processor.  Physically adjacent subdomains will need to communicate;
depending on the parallel architecture, connections may need to be
explicitly routed between their respective processors.  These
placement and routing steps are very important, because they determine
the efficiency of the runtime communication.  To minimize
communication time, adjacent subdomains should be placed on processors
which are separated by a short communication distance, and connections
should be routed so as to minimize the congestion of the most
congested physical links.  These are only general guidelines; the
properties of good placement and routing depend upon the architecture
of the target machine.

\noindent{\bf Time integration of semidiscrete equations.} After the
mesh has been partitioned, the next step is to integrate the system of
ordinary differential equations resulting from spatial discretization.
Integration schemes fall into two basic classes: explicit and
implicit. Explicit methods are conditionally stable but require only
sparse matrix-vector products at each time step. The simplicity of
their parallel implementation makes them attractive for many
hyperbolic problems \cite{ortega85}. Implicit methods, on the other
hand, necessitate solution of a sparse linear system at each time
step, but offer the possibility of larger time steps stemming from
their unconditional stability.  All else being equal, the choice of
integration scheme depends on the relative time and space scale of the
physics. For many wave propagation problems, choosing a time step for
accuracy often satisfies the stability requirement. However, for
interface problems with large material contrasts, the stability
requirement is often onerous, and implicit methods are preferable.

%At each time step the system of equations may be solved  either
%directly or by iteration.
%For high-frequency wave propagation through homogeneous regions, where
%the mesh size $h$ is dictated by length of incoming waves, and not
%geometric complexity or gradient of the displacement field, stability
%analysis of explicit methods yields time steps larger than those
%needed for accuracy. In such cases, explicit methods are more
%efficient.

%For wave propagation through inhomogeneous regions, explicit methods
%are most efficient when the mesh can be graded so that $h_{min}$ is no
%smaller than required. In this case, the stability
%criterion is never activated; accuracy controls. When $h_{min}$ is
%lower than needed to resolve the wave propagation, the time step
%becomes artificially low, and explicit methods suffer from too many
%iterations. In this sense structured meshes do double damage: finer
%grids imply more work per time step in addition to calling for a
%larger number of time steps. Even with unstructured meshes, the time
%step may be lower than required by the physics, especially for large
%contrasts in material properties, because the mesh resolution cannot
%change quickly enough, and smaller-than-desired elements ensue near
%transitions. Furthermore, geometric variability may occur in some
%regions at a smaller scale than the mesh needed to capture propagating
%waves, again resulting in smaller time steps than needed for accuracy.
%A third situation in which implicit methods may be preferable is for
%liquefaction problems in which large strain gradients may develop,
%again requiring finer resolution to capture. In all these cases, an
%implicit solver allows a larger time step, as dictated by accuracy, at
%the expense of having to solve an indefinite system of equations with
%sparsity structure defined by the mesh.

\noindent{\bf Parallel solution of linear algebraic systems.}
Implicit methods introduce the need to solve a linear algebraic system
at each time step, a task which can be avoided with explicit methods
by diagonalizing the mass matrix. Which is preferable depends on the
efficiency of a larger number of easily-parallelized sparse
matrix-vector products relative to a smaller number of linear system
solves. Parallel algorithms for the solution of sparse linear systems
of equations have received considerable attention; for example see
Heath et al.\ \cite{heath91}, Ortega \cite{ortega88} and Dongarra et
al.\ \cite{dongarra91}.  For linear waves, the coefficient matrix of
the linear systems is constant over time, so a direct method is ideal:
the factorization is performed just once, and only triangular solves
are required at each time step. Unfortunately, direct methods do not
map well to highly parallel systems. Current trends for large problems
favor variants of the conjugate gradient family (see Ashby et al.\
\cite{ashby90} for a framework), which minimize some norm of the
algebraic error over a Krylov subspace, and require preconditioning to
be effective for large problems.

\noindent Domain decomposition methods---solution of subproblems
obtained by partitioning a domain into subdomains, coupled with an
iterative resolution of subdomain interaction---have emerged as a
particularly effective method for parallel solution of elliptic
equations (e.g.\ \cite{Keyes,gropp88}.  The resulting reduction in
condition number from $O(1/h^2)$ for the overall system to $O(1/h)$
\cite{mansfield90} indicates many fewer iterations when conjugate
gradients are used as the iterative solver for the interface system,
with no increase in communication requirements. The solution of each
subproblem on a processor is easily performed using standard
sequential (e.g.\ sparse Gauss elimination) methods. Although a rich
body of theory and algorithms has been established for elliptic
problems (see the proceedings of the series of conferences, published
by SIAM \cite{ddm1,ddm2,ddm3,ddm5}, little attention has been paid to
hyperbolic problems (see \cite{quarteroni90,farhat93} for initial
work). Basic difficulties include how to partition unstructured meshes
into subdomains, and how to precondition the interface coefficient
matrix.

%For elastic wave propagation problems, the coefficient matrices of
%these linear systems can be derived from an equivalent positive
%definite, self-adjoint elliptic operator, so the implicit
%time-integrator essentially ``sees'' a sequence of elliptic problems.

%such as GMRES, as well as multigrid methods \cite{brandt77}.
%Multigrid, which uses a hierarchy of grids to smooth out errors, is
%particularly attractive because of its optimal serial complexity for
%many problems---$O(n)$ for $n$ variables. However, its parallel
%performance is hampered by the coarse grids: as much time is spent in
%communication as in computation
%\cite{mcbryan91}. Conjugate gradient methods, . Domain decomposition
%preconditioning, which increases the work to communication ratio on
%each processor in exchange for a reduction in iterations, has emerged
%as a particularly effective device for extracting coarse-grained
%parallelism. How such methods fare relative to traditional iterative
%techniques on computers with advanced communication systems remains to
%be seen.


\noindent {\bf Parallel code generation.} 
In addition to the other issues discussed so far, we have to specify
the program for the parallel computer.  To ease the development of
programs by application specialists, we want to allow the
specifications of such programs to solve PDE problems in a high-level
language, and automatically generate efficient parallel code from
programs written in that language.  There are several reasons why
parallel applications should be automatically generated from a
high-level specification, rather than written by hand.  First, few
scientists or engineers want to invest the amount of time necessary to
learn how to program multicomputers.  Second, a high-level language
allows one to quickly change the code and experiment with different
linear solvers, preconditioners, domain decomposition methods, ODE
integrators, etc.  Third, most parallel finite element applications in
the literature are written to solve a specific problem (e.g.\
Navier-Stokes equations), and are not general enough to use for other
domains.  A general-purpose finite element code generator will ensure
that the effort spent writing optimized parallel linear solvers can be
easily reused for other problems.  Fourth, all applications can be
ported to new parallel systems by porting the code generator.  Thus
the problem of parallel code generation becomes an issue of language
design and compilation, two areas on which there is an extensive
background in computer science\cite{lampsonhints}.

\noindent {\bf Mesh adaptivity.}
In the simplest cases, we refine previously generated meshes by
applying {\em a posteriori} error estimates. See the collected papers
in \cite{adaptive89} for both work on static and dynamic adaptive
meshing on sequential machines. However, runtime adaptivity on
parallel machines is much more difficult.  Data structures become more
complex, and dynamic load balancing becomes necessary~\cite{saltz90}.

\noindent {\bf Parallel system design.} To obtain good performance for
unstructured PDE problems, a parallel system (and its base software)
must satisfy a number of difficult requirements. These include fast
communication between processors, high memory bandwidth, good
compilers for the individual processors, and realistic and effective
tools for debugging and performance monitoring.

\noindent The most difficult issue is the need for fast irregular
communication between processors.  A general-purpose message passing
system that supports communication between arbitrary pairs of
processors in the system provides flexibility, but usually at the cost
of a high startup overhead.  Such overhead can reduce the effective
bandwidth by an order of magnitude if the application sends only small
messages with 10s of words.  Long-term connections between processors
(either via dedicated wires, e.g.\  in a hardware reconfigurable
computer, or by setting appropriate switches in a software
configurable computer) can reduce the communication startup overhead
at the expense of adding state information to the processor.  One of
the challenges is to find the right balance between flexibility
(little state), speed (lots of state), and the concerns of the
operating/runtime system (checking state).  CMU has been attacking
this problem for a number of years in the design of successive
generations of parallel systems.


%\noindent The most difficult issue is the need for fast irregular 
%communication between processors. CMU has been attacking this problem
%for a number of years. With its industrial partners, CMU has designed
%and built two parallel systems: the Warp~\cite{annaratone87e} and its
%successor, the iWarp~\cite{borkar88,borkar90}. Each system broke new
%ground in the design of communication systems. The Warp pioneered the
%idea of {\em direct access} to the communication system (systolic
%computing~\cite{kung82a}). The iWarp extended this idea with the
%notion of {\em direct access to programmable connections}.

%\noindent iWarp provides highly integrated hardware support that allows 
%a user program to rapidly and dynamically establish connections
%between arbitrary processors, without operating system overhead.  The
%result is an extremely efficient and flexible communication system
%that we can use as a laboratory for studying issues in fast
%interprocessor communication. To date, we have learned that we can use
%iWarp's communication system to exploit compile-time information about
%static communication patterns~\cite{FSW,stricker92}, develop highly
%efficient general-purpose message systems~\cite{subhlok93}, create
%efficient software techniques for barrier
%synchronization~\cite{feldmann92}, and provide an additional source of
%operands for the floating point units~\cite{gross92a}.  iWarp's
%communication mechanism is ideal for the irregular communications that
%arise in unstructured PDE problems.

\section{Preliminary Results}

The events that led to our collaboration and ultimately to this
proposal are interesting. In late 1991, the engineers (Bielak and
Ghattas) began exploring the use of commercial parallel computers for
modeling earthquake-induced ground motion. They implemented an
explicit solver on the CM-2 using an unstructured mesh. The code ran
slowly because the processors spent almost all of their time
communicating with each other over the relatively slow global
communication router.  So they rewrote the code, this time creating a
block-structured mesh and painstakingly partitioning it over the
processors to try to minimize the use of the relatively slow router.
The fast nearest-neighbor communication system was used for the
regular communication patterns within each block, while the global
communication system was invoked to ``stitch'' the blocks together.
The resulting code exhibited an order-of-magnitude improvement in
MFLOPS, due to the substantial reduction in global communication.
Unfortunately, because the gradation of the block-structured mesh was
insufficiently flexible, smaller elements than desired were generated
in the rock surrounding the valley. These elements were critical with
respect to the stability condition; thus, the solver was forced to
take many more iterations than the physics required. Switching to an
implicit solver was problematic, for it entailed a complete rewrite of
the code. Parallel computers, rather than a tool, were becoming an
impediment to doing science!
% this last sentence is folksy, and was put in as a joke, but i kinda
% like it! 

Frustrated by the painstaking level of detail required to create 
reasonably efficient programs, they sought out the computer scientists
(Gross and O'Hallaron) for assistance. Thus was born a collaboration
aimed at developing tools that allow researchers to model scientific
problems, to investigate different solvers, and to visualize results,
all while extracting maximum performance from distributed memory
machines without involvement in low-level issues. For the past year we
% we need a better term than "low level issues" (hardware issues?
% communication issues?) 
have worked closely together to understand the issues and to produce a
working prototype tool for solving unstructured PDE problems on
parallel systems. The early experience with the CM-2 helped us
understand the close dependence between performance, problem
formulation, solution technique, and the underlying architecture. In
the process we gained a deeper appreciation for the value of
higher-level automatic tools and the pressing need for such tools on
parallel systems. Below, we describe the initial fruits of our
collaboration---two prototype tools---and the hardware for which they
are well-suited. 

\noindent {\bf Compiler Technology: Archimedes.}
We have built an automated system, named Archimedes, for solving PDEs on
unstructured meshes using distributed memory
supercomputers~\cite{cacm93}.  We have implemented Archimedes on both the
iWarp~\cite{borkar88,borkar90} and the CM-5~\cite{CM5}.

\noindent Figure \ref{archimedes} diagrams the structure of Archimedes.
A user provides two things: a description of the surface geometry of the
problem domain, and a high-level algorithm for performing finite element
calculations which solve some physical problem.  Archimedes includes a
mesh generator, a mesh partitioner, placement and routing heuristics,
and a code generator.  With the exception of the mesh generator, all
of these components can be used for two or three-dimensional problems.

\noindent Our two-dimensional mesh generation algorithm is due to Ruppert
\cite{Ruppert}, and is built on the Delaunay triangulation primitives of
Guibas and Stolfi \cite{Guibas}. Ruppert has proven that the algorithm
can discretize two-dimensional regions with straight-line boundaries
using triangles with angles no smaller than $20^\circ$ (except where
there are small input angles; these cannot be improved).  In practice,
the algorithm produces relatively few nodes, as Figure \ref{mesh}
illustrates; typically five times fewer nodes than competing algorithms
such as quadtree methods \cite{BEG}.  Our mesh generator can also refine
meshes (based on {\em a posteriori} error estimates), as in Figure
\ref{refine}.  Extending these algorithms to three dimensions
is one component of our proposed research. 

\noindent Our partitioner is based on an algorithm designed at CMU by
Miller et al.\ \cite{MTTV92}.  The algorithm produces a provably
asymptotically optimal separator; that is, for an $n$-node mesh in
$d$-dimensional space, it produces separators of size $O(n^\frac{d-1}{d})$.
(A separator is a set of vertices which, when removed, separates a mesh
into two pieces of balanced size.)
In practice, it produces partitions which are competitive with more
expensive methods such as spectral bisection~\cite{PSL90}; however, unlike
spectral methods and neural-net based methods \cite{byun92}, it runs in
$O(n)$ time.  Furthermore, its stability does not degrade as $n$ becomes
very large, as is the case with spectral partitioners.  A sample
partition is illustrated in Figure \ref{partition}.  Figure
\ref{quotient} illustrates the communication graph of this partition.

\noindent Fast and effective placement and routing for iWarp have been
developed by Feldmann et al.\ \cite{FSW}.  Figure \ref{routing} illustrates
the pathways needed to route the communication graph on iWarp, using two
phases of communication.  The partitioner and the placement code share
the trait that they use geometric information from the mesh to find
fast solutions to otherwise difficult problems.  Although we have used
two-dimensional illustrations for clarity, we can also partition and place
three-dimensional meshes.

\noindent The user writes an algorithm for constructing and solving the 
finite element system for a specific problem.  Archimedes includes a
code generator that, like the PARTI compiler of Saltz et al.\
\cite{sussman92}, generates parallel code for efficient processing of
unstructured meshes.  Unlike PARTI, however, Archimedes programs are
written in C with built-in primitive operations tailored for finite
element methods, guided by engineering needs.  Archimedes allows the
user to write code at the level of elements and matrices, and insulates
the programmer from low-level data structures by offering primitives for
matrix operations, Gaussian quadrature over elements, matrix assembly,
boundary conditions, graphical output, and other functions.
This allows a researcher
to quickly and easily experiment with different solution techniques in
a way not possible with low-level compilers such as PARTI.
Archimedes programs are machine-independent, and can be
written without knowledge of the parallel machine's underlying
communication mechanisms.  Programs are compiled into code for
the iWarp or the CM-5.

\noindent We have used Archimedes to duplicate the 3D seismic ground 
motion study we originally did on the CM-2. The results have taught us
a number of things. We have learned that it is indeed feasible to
construct automatic tools for solving unstructured PDE problems on
parallel systems. We have also learned that it is possible for
automatic tools to generate efficient parallel code. For example,
because of careful attention to communication overheads, in both
Archimedes and in the design of the iWarp communication system, the
seismic ground motion simulation produced by Archimedes for the iWarp
system spends only 20 percent of its time transferring data between the
processors; the remaining 80 percent of the time is spent doing useful work.
Finally, we have learned that creating a communication system this
effective is not easy. For
example, a large sparse solver generated by Archimedes for the CM-5
spends 80 percent of its time transferring data between processors,
because of the relatively low bisection bandwidth of the
CM-5, and because operating system overhead on message passing is too
large for the small messages that naturally arise from irregular
communication patterns. Our conclusion is that careful attention must be
given to design of communication systems in future parallel computers.

\noindent {\bf Compiler technology: Parallelizing Fortran compiler.}
We have developed a parallel Fortran (HPF-like) compiler for
distributed memory systems~\cite{subhlok93}.  The primary target of
this compiler is the iWarp system, although we have also ported it to
the Intel iPSC/860 and Paragon systems.  The compiler performs many
``backend'' communication optimizations; high-level issues like data
placement and distribution are left to the user (specified using HPF
directives), but the compiler handles all details of data management.

\noindent Much of the irregular communication support pioneered in our
Fortran compiler has been or will be applied to Archimedes.  In
particular, we have developed a new message passing model called {\em
direct deposit}~\cite{subhlok93} which allows message data to be
scattered into memory on the receiving processor without buffering.
In iWarp test programs, the direct deposit method provides
communication time improvements of 30-80\% over conventional message
passing. We expect similar gains for Archimedes programs.

\noindent {\bf Parallel system design.}
CMU has a long tradition of impacting the design of parallel
computers; with its industrial partners, our research group has
designed and built two parallel systems: the Warp~\cite{annaratone87e}
and its successor, the iWarp~\cite{borkar88,borkar90}. Each system
broke new ground in the design of communication systems. The Warp
pioneered the idea of {\em direct access} to the communication system
(systolic computing~\cite{kung82a}). The iWarp extended this idea with
the notion of {\em direct access to programmable connections}.  During
the last couple of years, we have used iWarp as a platform to
characterize application requirements for parallel systems, based in
part on example problems generated by Archimedes and our Fortran
compiler.

\noindent iWarp provides highly integrated hardware support that allows 
a user program to rapidly and dynamically establish connections
between arbitrary processors, without operating system overhead.  The
result is an extremely efficient and flexible communication system.
To date, we have learned that we can use iWarp's communication system
to exploit compile-time information about static communication
patterns~\cite{FSW,stricker92}, develop highly efficient
general-purpose message systems~\cite{subhlok93}, create efficient
software techniques for barrier synchronization~\cite{feldmann92}, and
provide an additional source of operands for the floating point
units~\cite{gross92a}.  iWarp's communication mechanism is ideal for
the irregular communications that arise in unstructured PDE problems.

\section{Research agenda}

%outline:
%
%> overview
%  > goal is to do LA.
%  > will require work along three fronts:
%
%> math/phy modeling
%  > build the model
%  > validate the results
%
%> numerical approximation
%  > use tools developed below to experiment with existing methods 
%
%> computation
%  > tool building: Archimedes
%  > tool building and algorithm design: mesh adaptivity
%  > tool building: integrating Archimedes techniques with Fortran
%  > new system designs

With recent progress in the characterization of large basins through
seismic and geotechnical exploration techniques and with the rapid
advances in high-performance computing, the time is right to mount a
concerted effort for predicting strong seismic ground motion in large
US basins by computer simulation. We propose specifically to predict
surface ground motion during large earthquakes in the Greater Los
Angeles Basin, based on three-dimensional models that take into
consideration source, propagation path, and local site effects.  This
will greatly improve our ability to establish appropriate guidelines
for prescribing the spatial distribution of the design ground motion,
thereby contributing directly to the mitigation of seismic hazard to
the general population in this seismically active region.

Achieving the goal of predicting strong seismic ground motion in the
Los Angeles basin will require a unified effort in three areas:
physical and mathematical modeling, numerical approximation, and
high-performance computing. The basic approach in all cases will be to
try new ideas on a small well-characterized valley in the CIS (Chusal
Valley), and then to incorporate the successful ideas into the full
Los Angeles Basin model. The outcome of the research will be: (1) A
characterization and improved understanding of the ground motion of
the Los Angeles basin to a degree that has not been possible until
now; (2) a set of algorithms and software tools for solving
unstructured PDE problems on parallel systems, applicable to a wide
class of problems beyond the seismic applications considered here; and
(3) a characterization of the computation and communication
requirements of such systems, with a set of recommendations for the
design of future parallel systems.

We propose a four-year research program to attain these goals.
Following is a summary of the research agenda:

\subsection{Physical and mathematical modeling}
The Greater Los Angeles Basin is a large composite sedimentary basin
(about 140km$\times$100km$\times$20km) that consists of four individual
basins: San Fernando, San Gabriel, San Bernardino, and Los Angeles
proper (Figure \ref{lacontour} shows the San Gabriel, Los Angeles, and
a portion of the San Fernando basins). The individual basins are quite
complex in geometry and structure, and are embedded in Pre-Tertiary
bedrock and filled with Tertiary and softer Quaternary sedimentary
layers. The geological structure of the sediments has been
investigated extensively using oil drilling data
\cite{yerkes65,davis89}, and the subsurface shear-wave velocities have
been measured directly over the basin \cite{fumal85}.  One of the main
reasons for choosing the Los Angeles Basin as the main object of our
proposed research, besides being the most highly populated seismic
region in the country, is that it is a well-characterized geologic
structure (including a varied fault system), amenable to mathematical
modeling.  A number of earthquakes have been recorded recently in this
region.  Records of actual ground motion will be used to check our
numerical predictions. Also, whereas portions of the basin can be
expected to respond nonlinearly for sufficiently large excitations,
this behavior can be modeled by using available constitutive
elastoplastic relations \cite{Zien90a,dakoulas92}.  The more extreme
case of soil liquefaction, generally present only in saturated loose
sands, is unlikely to occur in this basin. Attenuation properties of
the various soil and rock materials will be incorporated into the
models by appropriate anelastic attenuation factors, Q, for analyses
in the frequency domain, and by introducing viscoelastic material
behavior \cite{carcione88,moczo93} in time domain analyses.

As an initial check of our methodologies we plan to analyze the
seismic response of Chusal Valley. Its small size (about
700m$\times$400m$\times$60m), simple shape and simple geologic
formation, plus the fact that its response to a number of earthquakes
is well-documented \cite{king84,tucker84}, make this valley an
excellent test site. Using our three-dimensional numerical codes, we
will also be able to compare against observations and numerical
results from previous two-dimensional analyses \cite{bard85}.

The next step will be to develop a simple three-dimensional model of
the Los Angeles basin, by itself, using a simplified geometry and only
a few horizontal layers to characterize the stratigraphy
approximately. We plan to analyze this single basin by finite elements
in the time domain. Due to their high accuracy and applicability to
medium-size problems, we will also use boundary integral methods
\cite{hisada93}, including a recently developed
three-dimensional Gaussian beam-boundary element method in the
frequency domain \cite{kato93}. A comparison of results from the two methods
for various sources with those from actual earthquakes will serve to
validate both techniques, and, by allowing inelastic behavior in the
finite element method, to investigate its influence on the seismic
response. We will then be able to proceed with increased confidence,
in a third stage, to building a more realistic model, both with
regards to geometry and to material properties, first of the Los
Angeles Basin, and subsequently of the entire Greater Los
Angeles Basin, using only the finite element method.

In all cases involving the use of the finite element method, it will
be necessary to introduce fictitious absorbing boundaries to limit the
domain of computation. A number of approximate techniques that help
minimize the occurrence of spurious reflected waves, in both halfspace
and layered systems, have been developed in recent years 
\cite{engquist79,bamberger88,kallivokas91,higdon92}.
Appropriate methods for incorporating the seismic excitation into the
model if the source lies outside the computational domain are also
available \cite{bielak84,cremonini88}.

A detailed study of the spatial and temporal distribution of the
seismic response for a number of earthquakes along different faults
will help resolve important outstanding questions in strong motion
seismology and earthquake engineering in large basins.  For instance,
is there a need for seismic microzonation for the design of
constructed facilities?; under what conditions can one use
one-dimensional flat-layered models for predicting the seismic
response at a given site?; to what extent can the effects of site
conditions be isolated from those due to source and path?; under what
conditions, if any, will inelastic soil behavior affect
significantly seismic response?; and, lastly, will a basement or layers
with folds (humps) produce a highly unusual focusing or defocusing of
seismic waves, resulting in localized or diffuse motion, resembling
chaotic behavior on the free surface?

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\subsection{Scientific computation and parallel systems}

Another goal is to develop tools and environments for solving PDE
problems with large, unstructured, meshes, using our preliminary work
as a basis.  Our theme will be to emphasize functionality for
earthquake ground motion modeling, while maintaining generality. In
particular, special attention will be given to developing the
following capabilities:

\noindent {\bf Tool building: Archimedes.}
We will extend the functionality of Archimedes in the following ways:
%, maintaining its generality while
%emphasizing the functionality needed to model the LA basin:
(1) We will build a three-dimensional Delaunay-based mesh generator
and investigate heuristics for improving the quality of the meshes.
(2) To enable us to work with very large meshes, we will parallelize
as much of the Archimedes toolchain as possible.  The partitioning
algorithm will be almost trivial to parallelize; a more challenging
research goal will be to design a (possibly partition-based) practical
parallel mesh generation algorithm.  (3) We will design language
abstractions that can easily represent typical finite element
operations and insulate the programmer from the hardware and low-level
data structures, without constraining either the programmer or the
compiler more than necessary. This is a difficult balance, and will
require continuing close cooperation between the finite element
experts and the compiler writers.  Extensions include language support
for designing and evaluating domain decomposition methods,
preconditioners, and multigrid methods.  (4) Finally, we will continue
to investigate improvements to the code generator. For example, as we
have a partitioner for unstructured meshes, we can use it to generate
a nested dissection ordering~\cite{Khaira} on each subdomain for use
in domain decomposition methods~\cite{Keyes}.  This ordering can also
be used to improve cache performance.  As another example, if we
choose a vector machine such as the CM-5 as a research base, we will
use new high-performance scan-based vectorization techniques, designed
at CMU, which have produced the world's fastest sparse matrix-vector
multiply~\cite{Zagha}.

\noindent Initially, will test our ideas on the iWarp systems at Carnegie
Mellon. Later, we will port Archimedes to one of the new parallel
systems, such as the Intel Paragon and its successors, or the Cray
T3D\footnote{To continue its collaboration with CMU, Intel has offered
to provide us with a Paragon system. This system will be installed
probably before the award date and then will become available in the
same manner as the iWarp systems (at no cost to this research
project). The Pittsburgh Supercomputing Center expects installation
of a T3D sometime in 1993; we expect to be able to use this system in
the second half of 1994.  Our budget is based on the assumption that
access to machines at the PSC will---as it has in the past---not incur
any additional expenses.}.

\noindent {\bf Development of parallel implicit solvers.} For a class
of problems (including highly heterogeneous wave propagation),
implicit methods are more effective {\em on sequential computers} than
explicit methods. Using domain decomposition ideas, we will develop
implicit methods that are well suited to parallel architectures.  Our
partitioner provides us with a mechanism to make asymptotically
optimal cuts.

\noindent {\bf Evaluation of existing numerical methods/algorithms.}
The availability of an easy-to-use tool such as Archimedes will allow
us to evaluate the performance of different parallel numerical methods
and algorithms on unstructured meshes, a task which must be performed
experimentally. For example, we will study the interplay among problem
size, number of processors, hardware platform, problem geometry, etc.,
on the relative performance of implicit vs.\ explicit time integrators,
different preconditioners, different drivers, domain decomposed vs.\
full space iteration, etc.

\noindent {\bf Tool building and algorithm design:  mesh adaptivity.} 
Our existing capability for adaptive mesh refinement handles the
problem of dynamic adaptation as a sequence of static problems. While
constituting a reasonable approach to steady problems, this is less
than ideal for time-dependent problems. In anticipation of problems
requiring dynamic adaptivity, we will investigate how to adaptively
refine meshes at runtime. To do this, we must consider questions of
mesh geometry, criteria for refinement and ``derefinement'', load
balancing, and dynamic sparse matrix data structures.  Also, we must
find an elegant way to integrate these into Archimedes' language.

\noindent {\bf Tool building:  integrating Archimedes with parallel
Fortran.} Archimedes has irregularly distributed data structures and
FEM-specific language extensions.  Our Fortran compiler has only
regular data distributions, but is more flexible.  Ultimately, we will
combine their functionality into a single compiler, thereby enabling
the solution of more complex problems. For instance, solution of
inverse boundary value problems and design optimization problems by
nonlinear optimization alternates between large, sparse matrix
computations and smaller, dense matrix solutions \cite{ghattas93}. The
development of such a compiler increases the applicability of our work
to scientists and engineers interested in a wider class of design and
modeling problems.

\noindent {\bf Parallel system design.}
In the long term, we need to influence the way parallel computers are
designed and programmed.  Careful consideration must be given to the
design of hardware for irregular communication, and to ensuring that
application programs can access the communication system without
substantial operating system overhead. Such impact on future parallel
systems will require time (since the systems that can reflect the
outcome of this research will not be commercially available until
after the end of our research period), but our arrangements with
industrial partners place us in a unique position to leverage the
results of our research.

\noindent We will continue to  work closely with Intel's Supercomputer
Systems Division, a major vendor of parallel systems, on the design
and requirements for future-generation parallel systems.  Our focus
will be the design of a fast communication system and the elements of
the software tool chain (e.g., compilers, program generators, runtime
system) that are required to allow an application fast and direct
access to the communication system.  I.e., there are many aspects of a
fast communication system, including transmission time on the wire,
routing time, buffer management time, application activation time,
synchronization overhead, scheduler overhead, etc., and we want devise
a hardware/software structure that allows the application programs to
manage as many aspects of this communication system as possible.  Only
by handing control of communication to the application program (and
the compilers used to translate the application program), will we be
able to reduce the overheads and make fast communication for irregular
computations a reality.  A first step in this direction is to devise a
model of message passing that includes all these aspects and that
allows compiler writers and parallel system designers to make rational
tradeoffs.  We will develop the model with special consideration given
to the design of parallel systems such as the Intel Paragon and its
successors.  We will use a collection of unstructured PDE problems
written with Archimedes to validate the model.

\subsection{Summary}
%Synergy??

Our research plans regarding the issues of physical and mathematical
modeling and the issues of scientific computation and parallel systems
are closely coupled. The activities of one area support the other, and
vice versa. At the end of the research period, we expect to be able to
simulate the ground motion effects of earthquakes in the Greater LA
basin, and this capability can be used in many ways, as sketched in
this proposal.  The results of our investigation into scientific
computing and parallel systems will enable this capability but are
expected to have an impact far beyond this specific problem.  A tool
that allows the application specialist to solve PDE problems creates
opportunities in areas as diverse as manufacturing process simulation,
environmental contamination, and cardiovascular blood flow. The
lessons learned from this research will be available to other computer
architects (in the form of a characterization of this application with
regard to its communication and computation requirements).
Furthermore, if we are successful, our research will result in the
inclusion of appropriate support for irregular communication in all the
parallel systems implemented by our industrial partner.


\section{Personnel and educational activities}

The proposed research will be conducted primarily by Jacobo Bielak and
Omar Ghattas, of the Department of Civil Engineering at Carnegie
Mellon University (CMU); Keiiti Aki, of the Department of Geological
Sciences, University of Southern California (USC); Thomas Gross and
David O'Hallaron of the School of Computer Sciences (SCS) at CMU; and
Francisco Sanchez-Sesma of the Institute of Engineering of the
National University of Mexico (UNAM). Brief biographies of these
faculty appear in later pages. In addition, Richard C.  MacCamy,
Professor of Mathematics and a member of the Center for Nonlinear
Analysis at CMU, with whom the PI has collaborated closely since
1982, will participate in the nonlinear analysis aspects of the
proposed research at no cost to the sponsor. Professor MacCamy is an
authority on boundary integral equations and has conducted extensive
research on the behavior of inelastic materials, nonlinear flow
through porous media, and nonlinear wave propagation. In addition to
faculty participation we seek support for six PhD students and one
Postdoctoral fellow: three in CMU-Civil Engineering, one of whom will
be co-advised by Sanchez-Sesma, two in CMU-SCS, and two at USC. One
additional student will work with Sanchez-Sesma at UNAM on this
project, at no cost to NSF.

All the members of the team will be actively involved in planning and
executing the proposed research. The physical and mathematical
modeling aspects will be the responsibility, primarily, of Aki,
Bielak, and Sanchez-Sesma. Together they will select the specific
physical models to be considered in the investigation and will develop
the corresponding mathematical models, with the participation of the
other members of the team. Aki will provide guidance on seismological
aspects of the work. As explained in the main body of the text, most
models will be analyzed both by the finite element method and with
boundary elements. The development of the finite element methodology,
including the supporting problem solving environment, will be a
primary task of the CMU group, while the USC group will concentrate on
the boundary element techniques. Boundary element techniques will also
be pursued at UNAM and CMU; to fully achieve an integrated effort
Sanchez-Sesma will spend three months a year at CMU as a visiting
professor.  The numerical approximation aspects of the research at CMU
will be the primary responsibility of Ghattas, while Gross and
O'Hallaron will lead on the computer science aspects of the
investigation.  Informal interaction between the various investigators
will occur continually during the course of the project, with Bielak
serving as coordinator.  A more formal annual meeting will be held at
CMU during the time of Sanchez-Sesma's visit, with the participation
of Aki and the rest of the USC team, to review the current work and to
plan for the coming year. Also, Bielak will visit USC and the Southern
California Earthquake Center (SCEC) once a year to attend its annual
meeting and help integrate the work on this project with the overall
work of the Center.

Members of the team have known and interacted with each other over a
period of time. Bielak and Sanchez-Sesma were colleagues for some
time, while Bielak was on the faculty at the Institute of Engineering
at UNAM. Aki and Bielak have been interested in examining the
significance of inelastic soil behavior in sedimentary soils during
strong earthquakes. In Fall 1991 Bielak participated, at Aki's
invitation, in a workshop at the SCEC on this topic. The CMU group has
been working together for a year and a half. In addition to specific
research results mentioned earlier, this collaboration between
engineers and computer scientists has already resulted in increased
interaction between students in the two disciplines.  Computer science
students now attend engineering courses in finite element methods and
applications, and have become interested in problems in engineering,
and mechanics students attend seminars on high-performance computing.
As a result of this collaboration, a weekly seminar series on
unstructured mesh problems in computational science has taken place
over the last semester. One of our major educational objectives is to
use our collaboration as a springboard for organizing a program in
computational science and engineering at CMU, with participation from
scientists, engineers, and computer scientists. The university
administration and colleagues outside our immediate group have
expressed an interest in such program. Our first activity along these
lines is the creation of a university-wide course on the practice of
scientific computing, jointly taught by computer science and civil
engineering faculty. 

The SCS at CMU already has in place several activities that support the
proposed research (which will, in turn, be reinforced by the proposed
projects). SCS was one of the first schools to introduce an
undergraduate class on programming parallel computers, which has been
taught once a year for the last 4 years. The School also offers a
highly selective Honors program in Computer Science that requires
students to perform a research project and to write an Honors thesis.
The proposed research provides a good framework for additional students
in the Honors program for research in parallel systems. Also, the
proposed research provides material for the undergraduate class by
allowing the faculty to expose the students to real problems.

The PI, a Hispanic American, has been successful in attracting
graduate students of Hispanic background to CMU. At least six of of them
have joined the Department of Civil Engineering or the the Department
of Engineering and Public Policy, and three have worked directly under
his supervision. We plan to make every effort to recruit for the
project at least one Hispanic-American graduate student. In addition,
we have admitted a woman to the graduate computational mechanics
program who is expected to work on this project.

\section{Industrial and government relationships}

The proposed research leverages directly several multi-year
government and industrial research projects at Carnegie Mellon.  The
iWarp architecture was developed under (D)ARPA funding (with
substantial cost-sharing by Intel), and research in program
development tools and their impact on computer architecture continues
to be supported. These tools substantially ease the task of writing
programs for this novel architecture.  CMU operates three iWarp
systems (one with 1/2 Gigabyte of memory.  The NECTAR prototype
network (a 100 Mbit/sec local area network developed with funding by
(D)ARPA) is currently evolving into the Gigabit NECTAR network.
Currently, a high-speed high-bandwidth HiPPI interface for iWarp (
developed by Network Systems Corporation with (D)ARPA funding) is
being deployed, and we expect to have a 100 MByte/sec HiPPI-interface
available by the end of the summer for all iWarp systems, allowing a
connection between the iWarp systems in the School of Computer Science
and the Cray YMP-C90 and the CM-2 at the Pittsburgh Supercomputing
Center (and a Cray T3D when it arrives later this year).

Carnegie Mellon has a long track record of impacting industry through
technology transfer.  The Warp and iWarp projects are good examples,
moving the idea of systolic computing into companies where they find
their way into commercial products.  Intel and Carnegie Mellon have
agreed to continue a collaborate on applications, tools, and system
design for Paragon and its follow-on systems.  These systems attempt
to leverage the lessons learned from iWarp into mainstream Intel SSD
processors and products.  Unstructured mesh methods and irregular
communication are being recognized as important domains for future
parallel systems, and so this research can have a direct and immediate
impact on the design of such future systems.

The proposed research also leverages a five-year industrial grant
recently awarded to Carnegie Melon's Computational Mechanics
Laboratory, whose director and associate director are the PI and one
of the Co-PIs. The research, sponsored by the finite element software
development company Algor, focuses on computational methodologies for
coupled, nonlinear field problems. One component of the research to be
undertaken under this project is the development of variational
boundary integral equation methods for nonlinear problems. If the
research is successful, its results will be applied to the nonlinear
earthquake ground motion problem. Another component of the research
will extend the methods for inducing sparsity in boundary integral
coefficient matrices, developed by the PI, to three-dimensional
problems. Again, these techniques have direct bearing on the proposed
research, and will be developed at no expense to NSF. The boundary
integral methodology we develop can be integrated easily into
Archimedes: support for dense matrix computations is a natural
outgrowth of our HPF work, and the ``sparsifying'' ideas of
\cite{zeng92} result in spatially-block dense matrices, whose
structure is identical to that of matrices produced by domain
decomposition methods, and thus are ideally suited to distributed
memory multiprocessors and to Archimedes, and will be supported. 

The proposed research will be integrated with the activities of the
Southern California Earthquake Center, whose director is one of the
Co-PIs.


\newpage

\bibliography{prop93}

\newpage

{\bf Jacobo Bielak} is a Professor of Civil Engineering at Carnegie
Mellon University. His primary research interests lie in the areas of
earthquake engineering and computational mechanics; specifically, on
the analysis of dynamic soil-structure interaction, inelastic
structural response, fluid-structure interaction, and seismic ground
motion, and on the development of numerical and computational
techniques for the efficient solution of these problems.

He is currently working on an \mbox{NSF-sponsored} exploratory project
on the \mbox{three-dimensional} ground motion modeling in sedimentary
basins using
finite elements on massively parallel computers; on the coupling of
boundary integral and finite element methods in wave propagation
interface problems; and on the development of artificial boundaries
for limiting the region of computation in finite element applications
of wave propagation problems.

The section on soil-structure interaction contained in the National
Earthquake Hazard Reduction Program Recommended Provisions for the
Development of Seismic Regulations for New Buildings is based
primarily on his work. He developed for the first time accurate
mathematical models to simulate on the Cray YMP computer the nonlinear
three-dimensional behavior of pile foundations. Most recently, he has
developed a novel methodology for generating sparse boundary integral
formulations for very large scattering problems, which have been
unmanageable till now.

Professor Bielak received his PhD in Civil Engineering from the
California Institute of Technology in 1971, where he also served as a
Research Fellow. From 1972 to 1978 he was on the faculty of the
Institute of Engineering of the National University of Mexico as Head
of its Dynamics Group. He joined Carnegie Mellon in 1978. He has
served as Chairman of the Dynamics Committee of the Engineering
Mechanics Division of the American Society of Civil Engineering, and
as a member of the editorial board of the Journal of Engineering
Mechanics. 

\begin{quote}
\footnotesize

J. Bielak, R.C. MacCamy, D.S. MacGhee, and A. Barry, ``Unified
Symmetric BEM-FEM for Site Effects on Ground Motion---SH Waves,''
{\em Journal of Engineering Mechanics}, Vol. 117, 2265--2285, 1991.

\vspace*{.5\baselineskip}

L.F. Kallivokas, J. Bielak, and R.C. MacCamy, ``Symmetric Local
Absorbing Boundaries in Time and Space,'' {\em Journal of Engineering
Mechanics}, Vol. 117, 2027--2048, 1991.

\vspace*{.5\baselineskip}

A.M. Trochanis, J. Bielak, and P. Christiano, ``Three-Dimensional
Nonlinear Study of Piles,'' {\em Journal of Geotechnical Engineering},
Vol. 117, 429--447, 1991. 

\vspace*{.5\baselineskip}

A. Barry, J. Bielak, and R.C. MacCamy, ``On Absorbing Boundary
Conditions for Wave Propagation,'' {\em Journal of Computational
Physics}, Vol. 79, 449--468, 1988.

\vspace*{.5\baselineskip}

J. Bielak and P. Christiano, ``On the Effective Seismic Input for
Nonlinear Soil--Structure Interaction Systems,'' {\em Earthquake
Engineering and Structural Dynamics}, Vol. 12, 107--119, 1984.

\end{quote}

\newpage
{\bf Omar Ghattas} is an Assistant Professor in the Department of
Civil Engineering and is affiliated with the NSF-Engineering Research
Center in Engineering Design at Carnegie Mellon. He is interested in
computational mechanics and optimal design, with special interest in
the development of numerical techniques for the solution of
large-scale optimization and computational mechanics problems.

In addition to collaboration with Prof. Bielak on seismic ground
motion simulation, his recent work has focused on the development of
parallel algorithms for optimal design, with application to shape
optimization of aerodynamic surfaces and aerospace vehicles, and on
the solution of these problems on a data-parallel computer.  He is
also working on topology optimization of solids and structures, and
finite element methods for constrained variational problems. 

Dr. Ghattas received the Ph.D. degree in Civil Engineering from Duke
University in 1988, where he also served as a research associate, and
joined Carnegie Mellon in 1989.

\begin{quote}
\footnotesize

C.E. Orozco and O.N. Ghattas, ``A Sparse Approach to
Simultaneous Analysis and Design of Geometrically-Nonlinear
Structures'', {\em AIAA Journal}, Vol. 30, No. 3, 1992. 

\vspace*{.5\baselineskip}

O.N. Ghattas and S.K. Schrader, ``A Nonlinear Programming
Formulation for the Topological Structural Optimization Problem'',
{\em Structural Optimization}, 1992.

\vspace*{.5\baselineskip}

X. Li, J. Bielak, and O. Ghattas, ``Three-Dimensional Earthquake
Site Response on a CM-2'', {\em Proceedings of the Tenth
World Conference on Earthquake Engineering}, Madrid, Spain, July 1992.

\vspace*{.5\baselineskip}

O.N. Ghattas and C.E. Orozco, ``A Parallel Algorithm for Shape
Optimization of Nonlinear Systems'', {\em Proceedings of the Fourth
Air Force/NASA/AIAA Symposium on Multidisciplinary Analysis and
Optimization}, Cleveland, Ohio, September 1992. 

\vspace*{.5\baselineskip}

C.E. Orozco and O.N. Ghattas, ``Optimal Design of Aerodynamic
Shapes on Massively Parallel Multiprocessors'', {\em Computing Systems
in Engineering}; also, {\em Proceedings of the Symposium on
High Performance Computing for Flight Vehicles}, Washington, DC,
December 1992.

\end{quote}

\newpage

{\bf Thomas R. Gross} is an Associate Professor in the School of
Computer Science at Carnegie Mellon University.  His primary research
interests lie in the area of computer systems, specifically the
tradeoffs between hardware, software, and architecture for the
implementation of high-performance systems.  He is currently working
on software tools (compilers, parallel program generators, runtime
systems, libraries) to simplify the task of programming parallel
computers.

While at Carnegie Mellon, he has worked (together with H.T. Kung and
other CMU researchers) on the Warp machine, the first programmable
systolic array that was programmed in a high-level language.  A
follow-on project, the iWarp architecture developed jointly with Intel
Corp., led to the development of the first programmable systolic array
that also supports more traditional models of inter-processor
communication (e.g., message passing).  This iWarp microprocessor
forms the core of the Intel iWarp product line.  Together with
Associate Dean Allan Fisher, Prof. Gross developed an innovative
course teaching the programming of parallel computers to
undergraduates.  This course has been offered for several years at
Carnegie Mellon and has been adopted by other universities.

Prof. Gross received his Ph.D. in Electrical Engineering from Stanford
in 1983, where he also received a MS degree in Computer Science.  He
also received the Diplom-Informatiker degree from the Technical
University in Munich, Germany.  His doctoral dissertation described
compiler optimizations and scheduling techniques for the
interlock-free pipeline structure of the MIPS machine, a reduced
instruction set microprocessor developed at Stanford University.
Thomas Gross joined Carnegie Mellon in 1984.

\begin{quote}
\footnotesize
S. Borkar, R. Cohn, G. Cox, T. Gross, H.T. Kung,
M. Lam, M. Levine, B. Moore, W. Moore, C. Peterson, J. Susman, 
J. Sutton, J. Urbanski, and J. Webb, ``Supporting Systolic and Memory
Communication in iWarp'', Proc. 17th Annual Intl.\ Symposium on Computer
Architecture, May, 1990, Seattle, pp.\ 70--81.

\vspace*{.5\baselineskip}
R. Cohn, T. Gross, M. Lam, and P. S. Tseng, ``Architecture and
Compiler Tradeoffs for a Wide Instruction Word Microprocessor''.  In
Proc.\ Third Intl.\ Conf.\ on Architectural Support for Programming
Languages and Operating Systems (ASPLOS III), pp. 2--14.  ACM/IEEE,
Boston, April 1989.

\vspace*{.5\baselineskip}
T. R. Gross, J. L. Hennessy, S. A. Przybylski, and C. Rowen,
``Measurement and Evaluation of the MIPS Architecture and Processor''.
ACM Trans. on Computer Systems 6(3):229--258, August 1988.

\vspace*{.5\baselineskip}
M. Annaratone, E. Arnould, T. Gross, H. T. Kung, M. S. Lam, O.
Menzilcioglu, and J. Webb, ``The Warp Machine: Architecture,
Implementation and Performance''.  IEEE Trans.\ on Computers
C-36(12):1523--1538, Dec.\ 1987.

\vspace*{.5\baselineskip}
J. Hennessy and T. Gross, ``Postpass Code Optimization of Pipeline
Constraints''.  ACM Trans.\ on Programming Languages and Systems
5(3):422-448, July 1983.
\end{quote}


\newpage
{\bf David R. O'Hallaron} is a Research Computer Scientist in the
School of Computer Science at Carnegie Mellon University. His research
interests are tools, applications, and architectures for parallel
systems.

Dr. O'Hallaron is experienced in diverse application areas, including
sonar digital signal processing, extended Kalman filtering, finite
element modeling, gravity modeling, tomographic imaging, general
matrix computations, and fast Fourier transforms. While at Carnegie
Mellon he has worked on the development of the Archimedes compiler and
the Fx parallel Fortran compiler, developed a parallel program
generator, now an Intel product, for building signal processing
applications, and led the effort to implement a real-time version of
the Navy's workhorse passive airborne sonar algorithm on iWarp.

He is currently involved in the development of parallelizing
compilers, as well as joint application efforts, with government,
industry, and academic partners, in the areas of earthquake site
response modeling, synthetic aperture radar imaging, passive sonar
processing, computational chemistry, and magnetic resonance imaging.

Dr. O'Hallaron received the Ph.D. degree in Computer Science from the
University of Virginia, Charlottesville, VA, in 1986.  From 1986 to
1989 he was with the GE Research and Development Center, Schenectady,
NY, as a Computer Scientist. In 1989, he joined the faculty in the
School of Computer Science at Carnegie Mellon University.\\

\begin{quote}
\footnotesize

\vspace*{.5\baselineskip}
J.~Subhlok, J.~Stichnoth, D.~O'Hallaron, and T.~Gross.  ``Exploiting
task and data parallelism on a multicomputer''.  Proceedings of
the ACM SIGPLAN Symposium on Principles and Practice of Parallel
Programming, San Diego, CA, May 1993.

\vspace*{.5\baselineskip}
A.~Feldmann, T.~Gross, D.~O'Hallaron, and T.~Stricker.  ``Subset
barrier synchronization on a private memory parallel system''.
Proceedings of the 4th Annual ACM Symposium on Parallel Algorithms and
Architectures, pp.~209--218, San Diego, CA, June 1992.

\vspace*{.5\baselineskip}
D.~O'Hallaron, ``Uniform approach for solving some classical problems
on a linear array''. IEEE Transactions on Parallel and Distributed
Systems, vol.~2, no.~2, pp.~236--241, April 1991.

\vspace*{.5\baselineskip}
D.~Nicol and D.~O'Hallaron, ``Efficient algorithms for mapping
pipelined and parallel computations''. IEEE Transactions on Computers,
vol.~40, no.~3, pp.~295--306, March 1991.
\end{quote}


\newpage
\pagestyle{empty}

\begin{figure}[p]
\vspace{7in}
\caption {\label{xli}
(a) Model basin; (b) Synthetic horizontal displacements; (c)
Distribution of horizontal displacements throughout the basin in E--W
direction at an instant; (d) Amplitude of peak E--W displacement
within the basin over entire history.}
\end{figure}

\begin{figure}[p]
\vspace{3.0 in}
\caption{\label{chusal}
Chusal response to local weak and strong motion (after
Tucker and King, 1984).}
\end{figure}

\begin{figure}[p]
\begin{minipage}{3.2in}
\psadobepercent{80}{kanto.color.ps}
\end{minipage}
\begin{minipage}{3.2in}
\psadobepercent{80}{contour.color.ps}
\end{minipage}
\caption{\label{kanto}
Interface between valley and basement in Kanto Basin, with
contour plot.}
\vspace{1in}
\end{figure}

\begin{figure}[p]
\centerline{\psadobepercent{150}{la.ps}}
\caption{\label{lasection}
A cross-section of the Los Angeles Basin}
\end{figure}

\begin{figure}[p]
\centerline{\psadobepercent{150}{la.struc.ps}}
\caption{\label{struc}
Structured mesh}
\end{figure}

\begin{figure}[p]
\centerline{\psadobepercent{150}{la.unstruc.ps}}
\caption{\label{unstruc}
Unstructured mesh}
\end{figure}

\begin{figure}[p]
\vspace{3 in}
\caption{\label{LA}
The Greater Los Angeles Basin (after Graves and Clayton,
1992).} 
\end{figure}

\begin{figure}[p]
\centerline{\psadobepercent{80}{archimedes.ps}}
\caption{\label{archimedes} Structure of Archimedes.}
\end{figure}

\begin{figure}[p]
\begin{minipage}{3.5in}
\psadobepercent{80}{mesh.ps}
\caption{\label{mesh} Finite element mesh of metal plate. \newline
\hspace{2in}}
\end{minipage}
\begin{minipage}{3.5in}
\psadobepercent{80}{refine.ps}
\caption{\label{refine} Mesh refined after stress simulation. \newline
\hspace{2in}}
\end{minipage}
\newline
\begin{minipage}{3.5in}
\psadobepercent{80}{part.color.ps}
\caption{\label{partition} Partitioned metal plate for \newline
sixty-four processor machine.}
\end{minipage}
\begin{minipage}{3.5in}
\psadobepercent{80}{quotient.color.ps}
\caption{\label{quotient} Sixty-four node communication \newline
graph for partitioned metal plate.}
\end{minipage}
\end{figure}

\begin{figure}[p]
\centerline{\psadobepercent{50}{routing1.ps}
\psadobepercent{50}{routing2.ps}}
\caption{\label{routing} Interprocessor connections on a 64-processor iWarp
for the communication graph.  Two communication phases are required.}
\end{figure}

\begin{figure}[p]
\vspace{4in}
\caption{\label{chusalcontour}
Plan view of portion of Chusal Valley, with
contour lines describing depth to basement (after King and Tucker,
1984).} 
\end{figure}

\begin{figure}[p]
\vspace{7in}
\caption{\label{lacontour}
Plan view of portion of the Greater Los Angeles Basin, with
contour lines describing depth to basement (after Yerkes, et al.,
1965).}  
\end{figure}
\end{document}

