Simulating the earthquake response of a large basin is accomplished by numerically solving the partial differential equations (PDEs) of elastic wave propagation, i.e. Navier's equations of elastodynamics. A variety of numerical methods have been used for approximating the solution of these problems. While boundary element methods have been popular for moderately-sized linear models, the inhomogeneity, nonlinearity, and large scale of such basins as the Greater Los Angeles Basin preclude their use here. On the other hand, uniform grid domain methods such as finite differences become impractical for the very large problem sizes involved.
To see why uniform grids are impractical, consider the Los Angeles
Basin. For a shear wave velocity of 0.4 km/s and a frequency of 2 Hz, a
regular discretization of the elasticity operator would place grid
points 0.02 km apart to achieve second order accuracy. The region of
interest has dimensions 140 km 
100 km 
20 km; thus, a
regular discretization, governed by the softest layer, requires 35
billion grid points with three displacement components per grid
point. At least a terabyte of primary memory would be needed, and on
the order of 
operations would be required at each time
step. The stability condition associated with explicit time
integration of the semidiscrete equations of motion imposes a time
increment at least as small as 0.004s. Thus, a computer would have to
perform at a sustained teraflop per second for two days to simulate a
minute of shaking.
Instead, we use unstructured mesh finite element methods that tailor the mesh size to the local wavelength of propagating waves. The wavelength is given by the product of the shear wave velocity and the period of excitation. The shear wave velocity is a property of the soil type; for a basin such as Los Angeles, it varies from 0.22 km/s to 4.5 km/s throughout the basin. Since in three dimensions mesh density varies with the cube of shear wave velocity, and since the softest soils are concentrated in the top layers, this means that an unstructured mesh method may yield three orders of magnitude fewer equations than with structured grids. Modeling the Los Angeles Basin for values of earthquake period and wave velocity that are desirable for engineering purposes thus becomes practical on the largest of today's parallel supercomputers.
We favor finite element methods for their ability to efficiently
resolve multiscale phenomena, the ease with which they handle stress
boundary conditions, and their firm theoretical foundation. For
temporal approximation, we have studied both explicit and
preconditioned conjugate gradient-based implicit methods. For
hyperbolic problems, explicit methods become unstable if the time step
is greater than the time it takes an elastic wave to cross any
element--the Courant condition. Implicit methods, on the other hand,
are unconditionally stable, implying that larger time steps can be
taken. However, the very characteristic that makes them stable--the
fact that the solution at time 
requires information from
all nodes
at time t--renders them unattractive on distributed memory
computers, since this implies global information exchange. We have
found that our mesh generators give us such good control over mesh
resolution that the Courant condition for explicit methods is not
onerous. The result is that the more readily parallelizable explicit
methods perform better for elastic wave propagation problems. In this
paper we discuss only a single-step explicit time integration method.
While unstructured mesh methods for simulating wave propagation through heterogeneous media result in many fewer equations, they introduce a number of computational difficulties that must be overcome. First, mesh resolution must closely follow wavelength; too coarse a resolution will introduce error, too fine will result in unnecessary computation as well as excessively small time steps (when explicit integration methods are used). Second, element aspect ratios must remain small; large aspect ratios will eventually result in instability in the time integration scheme. Highly heterogeneous basins, in which wavelengths vary rapidly in space, introduce special difficulties when trying to follow the wavelength change without severely stretching the mesh. Third, unstructured mesh methods are not easy to program on parallel computers; their irregular data structures require nontrivial mappings onto parallel machines, and irregular communication patterns are generated. Thus, we have had to develop fast, robust computational geometry and mesh generation techniques for highly spatially-variable meshes as well as compilers and tools that simplify the programming of unstructured mesh methods on parallel systems.
In the remainder of this paper, we present numerical methods, algorithms, and implementations for modeling earthquake-induced ground motion modeling in highly heterogeneous basins, and give performance results on the Cray T3D. We also describe Archimedes, a toolset/compiler we have built for automated solution of unstructured mesh PDE problems on parallel distributed memory computers. For an alternative approach to parallel ground motion modeling on distributed memory machines, see the work of Olsen, Archuleta, and Matarese [17], which employs finite differences on regular grids. See also the references to prior finite difference modeling work on sequential machines contained therein. In addition to the finite element method we describe in this paper, there have been recent efforts to endow finite difference wave propagation methods with multi-resolution capabilities. See the work described in [14], which uses composites of regular grids to achieve variable resolution.