Daniel Peter

Forward and adjoint simulations of seismic wave propagation on emerging large-scale GPU architectures

Computational seismology is an area of wide sociological and economic impact, ranging from earthquake risk assessment to subsurface imaging and oil and gas exploration. At the core of these simulations is the modeling of wave propagation in a complex medium. Here we report on the extension of the high-order finite-element seismic wave simulation package SPECFEM3D to support the largest scale hybrid and homogeneous supercomputers. Starting from an existing highly tuned MPI code, we migrated to a CUDA version. In order to be of immediate impact to the science mission of computational seismologists, we had to port the entire production package, rather than just individual kernels. One of the challenges in parallelizing finite element codes is the potential for race conditions during the assembly phase. We therefore investigated different methods such as mesh coloring or atomic updates on the GPU. In order to achieve strong scaling, we needed to ensure good overlap of data motion at all levels, including internode and host-accelerator transfers. Finally we carefully tuned the GPU implementation. The new MPI/CUDA solver exhibits excellent scalability and achieves speedup on a node-to-node basis over the carefully tuned equivalent multi-core MPI solver. To demonstrate the performance of both the forward and adjoint functionality, we present two case studies run on the Cray XE6 CPU and Cray XK6 GPU architectures up to 896 nodes: (1) focusing on most commonly used forward simulations, we simulate seismic wave propagation generated by earthquakes in Turkey, and (2) testing the most complex seismic inversion type of the package, we use ambient seismic noise to image 3-D crust and mantle structure beneath western Europe.

Anelastic sensitivity kernels with parsimonious storage for adjoint tomography and full waveform inversion

We introduce a technique to compute exact anelastic sensitivity kernels in the time domain using parsimonious disk storage. The method is based on a reordering of the time loop of time-domain forward/adjoint wave propagation solvers combined with the use of a memory buffer. It avoids instabilities that occur when time-reversing dissipative wave propagation simulations. The total number of required time steps is unchanged compared to usual acoustic or elastic approaches. The cost is reduced by a factor of 4/3 compared to the case in which anelasticity is partially accounted for by accommodating the effects of physical dispersion. We validate our technique by performing a test in which we compare the

A 1.8 trillion degrees-of-freedom, 1.24 petaflops global seismic wave simulation on the K computer

K_\alpha

Load-Balanced Local Time Stepping for Large-Scale Wave Propagation

sensitivity kernel to the exact kernel obtained by saving the entire forward calculation. This benchmark confirms that our approach is also exact. We illustrate the importance of including full attenuation in the calculation of sensitivity kernels by showing significant differences with physical-dispersion-only kernels.

Seismic Sound Lab: Sights, Sounds and Perception of the Earth as an Acoustic Space

We present high-performance simulations of global seismic wave propagation with an unprecedented accuracy of 1.2 s seismic period for a realistic three-dimensional Earth model using the spectral element method on the K computer. Our seismic simulations use a total of 665.2 billion grid points and resolve 1.8 trillion degrees of freedom. To realize these large-scale computations, we optimize a widely used community software code to efficiently address all hardware parallelization, especially thread-level parallelization to solve the bottleneck of memory usage for coarse-grained parallelization. The new code exhibits excellent strong scaling for the time stepping loop, that is, parallel efficiency on 82,134 nodes relative to 36,504 nodes is 99.54%. Sustained performance of these computations on the K computer is 1.24 petaflops, which is 11.84% of its peak performance. The obtained seismograms with an accuracy of 1.2 s for the entire globe should help us to better understand rupture mechanisms of devastating earthquakes.

Pillars of the Mantle: Imaging the Interior of the Earth with Adjoint Tomography

In complex acoustic or elastic media, finite element meshes often require regions of refinement to honour external or internal topography, or small-scale features. These localized smaller elements create a bottleneck for explicit time-stepping schemes due to the Courant-Friedrichs-Lewy stability condition. Recently developed local time stepping (LTS) algorithms reduce the impact of these small elements by locally adapting the time-step size to the size of the element. The recursive, multi-level nature of our LTS scheme introduces an additional challenge, as standard partitioning schemes create a strong load imbalance across processors. We examine the use of multi-constraint graph and hypergraph partitioning tools to achieve effective, load-balanced parallelization. We implement LTS-Newmark in the seismology code SPECFEM3D and compare performance and scalability between different partitioning tools on CPU and GPU clusters using examples from computational seismology.

Exascale Scientific Applications: Scalability and Performance Portability

We construct a representation of earthquakes and global seismic waves through sound and animated images. The seismic wave field is the ensemble of elastic waves that propagate through the planet after an earthquake, emanating from the rupture on the fault. The sounds are made by time compression (i.e. speeding up) of seismic data with minimal additional processing. The animated images are renderings of numerical simulations of seismic wave propagation in the globe. Synchronized sounds and images reveal complex patterns and illustrate numerous aspects of the seismic wave field. These movies represent phenomena occurring far from the time and length scales normally accessible to us, creating a profound experience for the observer. The multi-sensory perception of these complex phenomena may also bring new insights to researchers.

Variance-based Salt Body Reconstruction

In this work, we investigate global seismic tomographic models obtained by spectral-element simulations of seismic wave propagation and adjoint methods. Global crustal and mantle models are obtained based on an iterative conjugate-gradient type of optimization scheme. Forward and adjoint seismic wave propagation simulations, which result in synthetic seismic data to make measurements and data sensitivity kernels to compute gradient for model updates, respectively, are performed by the SPECFEM3D_GLOBE package [1] [2] at the Oak Ridge Leadership Computing Facility (OLCF) to study the structure of the Earth at unprecedented levels. Using advances in solver techniques that run on the GPUs on Titan at the OLCF, scientists are able to perform large-scale seismic inverse modeling and imaging. Using seismic data from global and regional networks from global CMT earthquakes, scientists are using SPECFEM3D_GLOBE to understand the structure of the mantle layer of the Earth. Visualization of the generated data sets provide an effective way to understand the computed wave perturbations which define the structure of mantle in the Earth.

Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes

Author(s): Straatsma, TP; Antypas, KB; Williams, TJ | Abstract:

Structure of the European upper mantle revealed by adjoint tomography

Seismic inversions of salt bodies are challenging when updating velocity models based on Born approximation- inspired gradient methods. We propose a variance-based method for velocity model reconstruction in regions complicated by massive salt bodies. The novel idea lies in retrieving useful information from simultaneous updates corresponding to different single frequencies. Instead of the commonly used averaging of single-iteration monofrequency gradients, our algorithm iteratively reconstructs salt bodies in an outer loop based on updates from a set of multiple frequencies after a few iterations of full-waveform inversion. The variance among these updates is used to identify areas where considerable cycle-skipping occurs. In such areas, we update velocities by interpolating maximum velocities within a certain region. The result of several recursive interpolations is later used as a new starting model to improve results of conventional full-waveform inversion. An application on part of the BP 2004 model highlights the evolution of the proposed approach and demonstrates its effectiveness.