Solving 3D Anisotropic Elastic Wave Equations on Parallel GPU Devices |
Any numerical implementation of equations 1-3 requires specifying both a computational mesh and a numerical discritization scheme. We define a computational grid and use time steps such that a grid point in space and time can be represented by quadruplet , where the integer counters have the ranges , , and . Continuous wavefield displacements at a given location are represented on the discritized grid as
(4) |
(6) |
starOfCubes
Figure 1. The 25 data points required to approximate a first derivative at the center shaded point using a FD stencil of 8 -order accuracy. |
---|
These difference operators allow us to specify a time-stepping scheme to calculate wavefield displacements in 3D elastic anisotropic media throughout the model domain. However, additional work is required to treat both the free-surface boundary and to minimize the energy in non-physical exterior boundary reflections. Implementing free-surface boundary conditions is fairly straightforward on a uniform grid because a (topography-free) free surface can be placed directly on the boundary which, assuming the free-surface normal vector points in the -direction, allows us to set where . We treat the other boundaries using a cascade of two operators comprised of absorbing boundary conditions (ABCs) derived from a one-way wave equation (Clayton and Enquist, 1977) and an exponential-damping sponge layer (Cerjan et al., 1985) of at least 48 grid points.
Figure 2 presents the nine procedural steps in the 3D FDTD algorithm required for calculating each forward time-step: 1) Compute strains from wavefield displacements according to equation 1; 2) Calculate stresses from constitutive relationship and according to equation 2; 3) Enforce the free-surface boundary condition (if required); 4) Inject a stress source (if not an acceleration source); 5) Compute acceleration from stress tensor (i.e., ); 6) Inject an acceleration source (if not a stress source; ); 7) Compute the forward time step from current and previous wavefield values; 8) Apply boundary conditions through cascade of ABC and sponge operators; and 9) Output data/wavefield as required. The next section details our GPU implementation of these steps.
flow
Figure 2. Control flow for FDTD algorithm implemented on a single GPU, including the division of tasks between the CPU host and GPU device. |
---|
Solving 3D Anisotropic Elastic Wave Equations on Parallel GPU Devices |