In this section we provide timings that characterize the performance of our parallel explicit wave propagation code on the Cray T3D. We are currently using the code to study the earthquake-induced dynamics of the San Fernando Valley in Southern California. The San Fernando simulations involve meshes of up to 77 million tetrahedra and 40 million equations. The largest mesh corresponds to the case of a lowest shear wave velocity of 220 m/s and a maximum frequency of 2 Hz; the code requires nearly 16 Gb of memory and takes 6 hours to execute for 40,000 times steps on 256 processors of the Cray T3D at the Pittsburgh Supercomputing Center (PSC). Results of a typical simulation, in which the basin was subjected to a vertically-incident plane wave Ricker pulse with a central frequency of 1 Hz, indicate a factor of amplification of eight in the softer parts of the basin compared to the maximum displacement on rock. This suggests substantially greater damage in these regions. A typical result is shown in Figure 7, which depicts the amplification induced by the soft soil. Simulations of this type are essential to better predict the local site effects within soft basins such as those on which Los Angeles, San Francisco, Mexico City, and Tokyo are situated.
Figure 7:
Surface distribution of ground motion amplification factors
in the San Fernando Valley. The amplification factors have been
calculated by comparing the surface to the bedrock motion.
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 our earthquake ground motion simulation 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. Additionally, the mesh sf1b corresponds to a geological model that 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 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 1.
Table 1: 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 input the file produced by Parcel, as well as the time to output results every tenth 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 8 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.
Figure 8: 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.
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 Figure 8.
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
Equation 6 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.
Let us extend this analysis to explicit methods on unstructured
meshes. Suppose 
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 of
processors. How close do we come to this ideal situation? First, we
plot the log of the computation to communication ratio against the log
of the number of processors, using the data in Figure 8. A
least-squares fit yields a line with slope 
. For a regular grid
with perfect partitioners, we have seen in the previous paragraph that
this slope should be 
. This suggests that the idealized
analysis is roughly applicable here.
Next, 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
9 shows the aggregate performance of our code
on the T3D
Figure 9: 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.
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 heterogeneous 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.
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 Equation
6--vector sums, diagonal matrix inversions, and sparse
matrix-vector multiplies require 
operations.
Figure 10 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.
Figure: 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 Equation
5 so that only one matrix-vector product is
performed at each time step.
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.