Whereas the mesh generation, partitioning, and parceling steps are currently performed sequentially, the wave propagation equations are solved on a parallel machine. This section describes the numerical techniques we use. Navier's equations of elastodynamics for an isotropic, heterogeneous medium are
where u is the displacement vector field, is the density, and and are elastic material constants, which depend on the shear and dilatational wave velocities.
The domain of the problem is an elastic halfspace, i.e.\ semi-infinite. In order to render the computational domain finite, we introduce absorbing boundaries at the bottom and sides of the basin that are local in both space and time [4, 10]. These boundaries allow the passage of outgoing waves with minimum reflection.
Since, in many cases, the earthquake source can be outside the computational domain, its effect must be introduced into the region. This is carried out as described in [3, 6] by means of effective forces. In short, for an arbitrary earthquake excitation these forces are determined in terms of the free-field motion by introducing a fictitious auxiliary surface that surrounds the basin. Across this auxiliary surface one imposes the conditions of continuity of displacement and traction. By selecting the total displacement vector field as the unknown in the resulting interior region and the scattered displacement field in the exterior region, the free-field displacement and traction now appear explicitly in the continuity conditions, which become jump conditions, with the free-field displacement and traction on the righthand side. These non-homogeneous terms on the righthand side are the ones that give rise to the effective forces upon spatial discretization. If, on the other hand, the seismic source is located inside the computational domain, say as a kinematic dislocation across the fault, one can select the fault itself as the auxiliary surface. The procedure is similar, but now one uses the total displacement everywhere as the unknown field; thus, the displacement field again experiences a jump across the interface, but the traction remains continuous. Notice that with this technique, whether the source is originally located inside or outside the computational domain, only outgoing waves will impinge upon the absorbing boundary. Both types of source are implemented in our code.
We also model material damping in the basin via viscous damping. With these modifications, standard Galerkin discretization in space by finite elements produces a system of ordinary differential equations of the form
where M is the mass matrix, C is the damping matrix associated with the absorbing boundary and material damping, K is the stiffness matrix, and f is the effective force vector. Here, M, C, and K are block matrices; the th block of M is a matrix given by
and the th block of K is given by
where is the finite element global basis function associated with the i-th node.
Damping is introduced through a proportional damping approximation at the element level, i.e. we take
where and are scalar constants and the superscript e indicates an element matrix. The first term leads to as damping factor that is inversely proportional to frequency, and the second to one that is linear in frequency. The constants and , which may vary within the basin according to the type of material, are chosen to best fit a prescribed attenuation law.
Given appropriate initial conditions, this system of ODEs can be integrated in time using central differences, yielding the explicit method
This method exhibits second-order accuracy in time, and when coupled with linear finite elements, we obtain second-order accuracy in space as well. We use a lumped mass approximation to M, which amounts to numerically integrating (3) with integration points at element vertices. This results in a diagonal mass matrix. To render the lefthand side operator of (6) diagonal, we further evaluate the off-diagonal components of at rather than . Inversion of the time stepping operator thus requires only a scaling of the righthand side of (6), which is carried out just once prior to time stepping. Forming the products of K and C with vectors comprises the major computational effort associated with iterating on (6). The sparsity structure of K is dictated by the underlying finite element mesh, and is thus very irregular. If shear waves are not over-resolved, the time step necessitated by stability is of the order of the time step dictated by accuracy, which is what an implicit method would take. By choosing an explicit method we avoid solving linear systems at each time step. Thus, overall, the explicit method is superior for our application, especially on a parallel computer.
We have tried several different choices of basis function order, and have concluded that piecewise-linear functions are the most efficient for problems requiring engineering accuracy. Our conclusion is based on numerical experimentation using plane Ricker wavelets on unstructured homogeneous meshes (in which case we know what the exact solution should be), but a simple argument can be given as follows. We recognize first that (spatial) approximation errors are bounded from above by interpolation errors. We then ask, for a given order of basis function and a given acceptable level of infinity-norm error, how many nodal points are required to produce a piecewise-polynomial interpolant of a simple harmonic wave. Next, we convert the required number of nodes per wavelength to an estimate of the storage and work required for an iteration of the explicit method (6). For example, if N is the total number of nodes, one can show that trilinear hexahedra require words of storage and 498N flops/time step, while triquadratic hexahedra necessitate 367N words and 1164N flops/time step. So, for example, triquadratic elements should require at least 2.2 times fewer nodes in order to be preferred (for storage reasons) over trilinear elements. However, one-dimensional interpolation theory tells us that 5% error requires 10 nodes per wavelength using linear elements or 9.4 nodes using quadratics. Thus, in three dimensions, triquadratics only allow times fewer nodes than trilinears, and are thus not warranted. An opposite conclusion is reached if one demands 99% accuracy. Our confidence in the values of material properties and in the fidelity of the source models for this problem does not warrant solution accuracies greater than 95%. Thus, we conclude that for this level of accuracy, the powers of higher-order interpolation are offset by their increased cost, both in storage and in increased work.