Efficiently modeling seismic datasets in complex 3D anisotropic media by solving the 3D elastic wave equation is an important challenge in computational geophysics. Using a stress-stiffness formulation on a regular grid, we present a 3D finite-difference time-domain (FDTD) solver using a 2
-order temporal and 8
-order spatial accuracy stencil that leverages the massively parallel architecture of graphics processing units (GPUs) to accelerate the computation of key kernels. The relatively small memory of an individual GPU limits the model domain sizes that can be computed on a single device. To circumvent this constraint and move toward modeling industry-sized 3D anisotropic elastic data sets, we parallelize computation across multiple GPU devices by using domain decomposition and, for each time step, employing an inter-device communication protocol to exchange data values falling within interior boundaries of each sub-domain. For two or more GPU devices within a single compute node, we use direct peer-to-peer (i.e., GPU-to-GPU) communication, while for networked nodes we employ message-passing interface (MPI) directives to route data over the network. Our 2D GPU-based anisotropic elastic modeling tests achieved a 10
speedup relative to an OpenMP CPU implementation run on an eight-core machine, while our 3D tests using dual GPU devices produced up to a 28
speedup. The performance boost afforded by the GPU architecture allows us to model seismic data for 3D anisotropic elastic models at lower hardware cost and in less time than previously possible.