3D modeling of time-harmonic seismic waves via a massively parallel structured multifrontal solver and applications
We consider the discretization and approximate solutions of equations describing time-harmonic seismic waves in 3D. We discuss scalar qP polarized waves, and multi-component elastic waves in inhomogeneous anisotropic media. The anisotropy comprises general (tilted) TI and orthorhombic symmetries. Particularly we are concerned with solving these equations on a large domain and for a large number of forcing terms, in the context of seismic problems in general, and modeling in particular. We consider variable order finite difference schemes, to accommodate anisotropy on one hand and allow higher order of accuracy on the other. ^ We resort to a parsimonious mixed grid finite difference scheme equipped with Perfect Matched Layer (PML) boundary conditions, for discretizing the Helmholtz operator, resulting in a non-Hermitian, indefinite and ill-conditioned matrix. We make use of a nested dissection based domain decomposition with separators of variable thickness, and introduce an approximate structured multifrontal solver by developing a parallel Hierarchically SemiSeparable (HSS) matrix compression, ULV factorization, and solution approach. The HSS tree based parallelism is fully exploited at the coarse level. The BLACS, PBLAS, and ScaLAPACK libraries are used to facilitate the parallel dense kernel operations at the fine-grained level. The communication patterns are composed of intra-context and inter-context ones. We exploit a randomized sampling approach to reshape the HSS construction phase, resulting in the extend-add of tall skinny sampling matrices, rather than the one associated with dense update matrices. The computational complexity associated with the entire structured multifrontal factorization is almost linear in the size, n say, of the matrix, viz. between O(nlog n) and O( n4/3 log n), while the storage is almost linear as well, between O(n) and O(nlog n). ^ We present various examples demonstrating the performance of our algorithm, as well as applications in multi-component modeling, inverse scattering transform or the so-called artifact free Reverse Time Migration (RTM), and illumination analysis. We present a comprehensive framework for wave-equation illumination analysis and introduce a target-oriented illumination correction that simultaneously accounts for limited acquisition aperture, and locally compensates for the so-called normal operator in inverse scattering to yield a partial reconstruction of reflectivity or reflection coefficient, while minimizing orientation dependent phase distortions and artifacts. To carry out the analysis we make use of higher-dimensional curvelets, which provide the means of extracting directional information.^
Maarten V. de Hoop, Purdue University.