We next discuss the performance of our parallel explicit wave propagation code on the Cray T3D. The relevant scenario for assessing the performance of our earthquake simulations as the number of processors increases is one in which the problem size increases proportionally, because unstructured PDE problems are typically memory-bound rather than compute-bound. Given a certain number of processors, we typically aim at full use of their memory; as the number of processors increases, we take advantage of their additional memory by increasing the problem size. In order to study the performance of the code with increasing problem size, we have generated a sequence of increasingly-finer meshes for the San Fernando Basin. These meshes are labeled sf10, sf5, sf2, and sf1, and correspond to earthquake excitation periods of 10, 5, 2, and 1 second, respectively, and a minimum shear wave velocity of 500m/s. Additionally, the mesh sf1b corresponds to the geological model used for the simulations described in the preceding section, which includes much softer soil in the top 30m, and thus necessitates an even finer mesh. Note that mesh resolution varies with the inverse cube of the excitation period, so that halving the period results in a factor of eight increase in the number of nodes. Characteristics of the five meshes are given in Table 2.
Table 2: Characteristics of San Fernando Basin meshes.
Our timings include computation and communication but exclude I/O. We exclude I/O time because in our current implementation it is serial and unoptimized, and because the T3D has a slow I/O system. I/O time involves the time at the beginning of the program to load and read the input file produced by the parceling operation, as well as the time to output results every 30th time step to disk. With the availability of the Cray T3E at PSC, we plan to address parallel I/O in the future.
We begin with a traditional speedup histogram, for which the problem
size is fixed and the number of processors is increased.
Figure 13 shows the total time, as well as the relative
time spent for communication and computation, for an earthquake ground
motion simulation, as a function of the number of processors. The
mesh used for these timings is sf2. On 16 processors, the time
spent for communication is 5% of the time spent for computation,
which is quite good for such a highly irregular problem. There are
about 24,000 nodes per processor, which results in about half the
memory on each processor being used. As the number of processors
doubles, the percentage of time spent communicating relative to
computing increases, as expected. For 128 processors, the
communication time has increased to one-fifth of the total time.
However, we are only utilizing 
of the local memory on a
processor; practical simulations will generally exhibit performance
more like the left bar of Fig. 13.
Figure 13: Timings in seconds
on a Cray T3D as a function of number of processors (PEs), excluding
I/O. The breakdown of computation and communication is shown. The mesh
is sf2, and 6000 time steps are carried out.
We can quantify the decrease in computation to communication ratio for
a regular 
mesh. Suppose there
are 
nodes on a processor, where P is the number of
processors. Suppose further that the regular grid is partitioned into
cubic subdomains of equal size, one to a processor. Since computation
for an explicit method such as Eq. 7 is proportional to
the volume of nodes in a cube (subdomain), and communication is
proportional to the number of nodes on the surface of the cube, the
computation to communication ratio is proportional to
, i.e. the ratio of total nodes to surface
nodes of the cube. Thus, for fixed N, the ratio is inversely
proportional to 
, at least for cubically-partitioned regular
grids with large enough numbers of nodes per processor. Clearly, it is
in our interest to keep 
as large as possible, if we want to
minimize communication time.
Consider now the case of unstructured, rather than regular,
meshes. Suppose that 
remains constant for increasing N and
P, i.e. the number of nodes per processor remains constant. Now
suppose that we have a partitioner that guarantees that the number of
interface nodes remains roughly constant as N and P increase
proportionally. Then we can expect that the computation to
communication ratio will remain constant as the problem size
increases. In this case, we have a method that scales linearly: the amount of time
required to solve a problem that is doubled in size is unchanged if we
double the number OCof processors. Let us attempt to hold the number of nodes per processor roughly
constant, and examine the aggregate performance of the machine as the
problem size increases. It is difficult to maintain a constant value
of 
, since processors are available in powers of two on the T3D.
However, we can still draw conclusions about scalability.
Figure 14 shows the aggregate performance of our
code on the T3D in megaflops per second, as a function of number of
processors (and, implicitly, problem size). Megaflops are those that
are sustained by matrix-vector product operations (which account for
80% of the total running time and all of the communication) during a
San Fernando simulation, exclusive of I/O. This figure shows nearly
ideal scalability, which is defined as the single processor
performance multiplied by the number of processors. These results show
that excellent performance is achievable, despite the highly
multiscale mesh. This behavior requires a partitioner that keeps the
number of interface nodes relatively constant for problem size that
increases concomitantly with number of processors.
Figure 14: Aggregate performance on Cray T3D
as a function of number of processors (PEs). Rate measured for
matrix-vector (MV) product operations (which account for 80% of the
total running time and all of the communication) during 6000 times
steps.
An even better measure of scalability is to chart the time taken per
time step per node. If the algorithm/implementation/hardware
combination is scalable, we expect that the time taken will not change
with increasing problem size. Not only must the partitioner produce
``scalable'' partitions for this to happen, but in addition the PDE
solver must scale linearly with N. This happens when the work per
time step is 
. This is obvious from the iteration of
Eq. 7--vector sums, diagonal matrix inversions, and
sparse matrix-vector multiplies require 
operations.
Figure 15 depicts the trend in unit wall clock time as the number of processors is increased. Unit wall clock time is measured as microseconds per time step per average number of node per processor, which includes all computations and communications for all time steps, but excludes disk I/O. As we have said above, for a truly scalable algorithm/implementation/hardware system, this number should remain constant as problem size increases with increasing processors. The figure demonstrates that we are close to this ideal. Ultimately, wall clock time per node per time step is the most meaningful measure of scalable performance for our application, since it is a direct indicator of the ability to solve our ultimate target problems, which are an order of magnitude larger than the San Fernando Basin problem we have described in this paper.
Figure 15: T3D wall-clock time in microseconds per time step per
average number of nodes per processor (PE), as a function of number of
processors. This figure is based on an entire 6000 time step
simulation, exclusive of I/O. The sf1b result is based on a
damping scheme in which 
in Eq. 6 so
that only one matrix-vector product is performed at each time step.
