Manuel Baumann

Fast Iterative Solution of the Time-Harmonic Elastic Wave Equation at Multiple Frequencies

Manuel M. Baumann

This work concerns the efficient numerical solution of the elastic wave equation. The elastic wave equation is a well-established partial differential equation (PDE) that models wave propagation through an elastic medium such as the earth subsurface, and is, therefore, of great importance in seismic applications. Geophysicists match simulation results of the elastic wave equation with measurements in a PDE-constrained optimization framework in order to gain information about the structure of the earth subsurface. After spatial discretization and Fourier transform in time, the time-harmonic elastic wave equation (forward problem) reads,

(K + iωkC − ωk2M )xk = b, k = 1, ..., Nω, (#)

where the primary challenge of this work is the efficient numerical solution of (#) when multiple (angular) wave frequencies are present, that is Nω > 1. This task becomes in particular challenging when the elastic wave equation in three spatial dimensions is considered because the matrices K, C and M in (#) become very large and ill-conditioned. In this situation, Krylov subspace methods are the common choice for the iterative numerical solution of (#). Without appropriate preconditioning, however, the Krylov iteration converges slowly to the solutions of the linear systems (#). The main contributions of this work are:

1. Development of an efficient shift-and-invert preconditioner designed for the simultaneous iterative solution when (#) is reformulated as a sequence of shifted linear systems. The shift-and-invert preconditioner is optimal with respect to a spectral convergence bound of multi-shift GMRES.

2. Implementation and development of an algorithmic framework that solves shifted linear systems in a nested inner-outer iteration loop. The new algorithm has been evaluated for different combinations of inner and outer multi-shift Krylov methods.

3. In a practical application, the preconditioner is usually applied inexactly. We, therefore, extend the recent theory of multilevel sequentially semiseperable (MSSS) matrix computations to the elastic operator in two and three spatial dimensions.

This work can be seen as a continuation of the extensive research of the last decade on the Complex Shifted Laplace preconditioner for the (discretized) acoustic wave equation that has been performed to a large extend at Delft University of Technology. The other way around, many of the contributions of this thesis also apply to the Helmholtz operator in a multi-frequency setting. We conclude this work with various numerical experiments in two (d = 2) and three (d = 3) spatial dimensions. In both cases, the size of the computational mesh in one spatial direction n is restricted by the highest wave frequency considered such that typically 20 points per wave length are guaranteed. If > 0 denotes the viscous damping parameter, our algorithm has shown computational complexity O(n^{d+1}) when the grid size and the frequency range are increased simultaneously, and multiple frequencies within this range are present.

Preconditioning techniques rely to a large extend on (numerical) approximations. The approximation can, in general, be motivated by physical insight of the dynamical behavior of the underlying PDE, or purely depend on the algebraic structure of the matrix that is obtained after discretization. In this work we have exploited both approaches to some extend: Inexact MSSS matrix computation techniques limit the growth of the off-diagonal sub-matrix rank and, hence, rely on the structure of the discretization matrix on a Cartesian grid. On the other hand, physical insight is used when the shift-and-invert preconditioner is solved at a damped frequency. Due to the damping, it is possible to efficiently slice a 3D problem into a sequence of 2D problems (block SSOR preconditioner) and additionally solve the problem on a smaller 3D grid (additive coarse grid correction).