Finite Element Methods for Seismic Imaging: Cost Reduction through mass matrix preconditioning by defect correction

Demand for hydrocarbon fuel is predicted to keep increasing in the coming decades in spite of easily accessible alternative fuels due to shifting geopolitical and economic situations. In order to find new hydrocarbon pockets, we need sharper images of earth’s subsurface. Also, the exploration of other sources of energy like geothermal will benefit from better models of the what lies underneath the surface. One way to obtain better images is to use superior numerical methods for forward modelling - Finite-element methods (FEM) are one such method, but their accuracy comes at the cost of increased compute expense. This thesis explores means to reduce this cost and adapt FEM to large-scale problems in geophysics. The Finite Difference (FD) method is the most popular numerical approximation scheme used in subsurface imaging problems. Representing the wave as the solution of individual motion and material equations is advantageous in terms of accuracy and stability and leads to the natural inclusion of density variations in the medium. This representation is referred to as the first-order formulation of the wave equation in this document. Finite Element (FE) methods are commonly derived for second-order equations because of the nature of variational formulation. Finite-element discretisations of the acoustic wave equation in the time domain often employ mass lumping to avoid the cost of inverting a large sparse mass matrix. Unfortunately, for a first-order system of equations, mass lumping destroys the superconvergence of numerical dispersion for odd degree polynomials. In chapter 3 of this thesis, we consider defect correction as a means to restore the convergence. We adapt the defect correction method to FEM by solving the consistent mass matrix with the lumped one as preconditioner. For the lowest-degree element, fourth-order accuracy in 1D can be obtained with just a single iteration of defect correction. In this chapter, we analyse the behaviour of the error in eigenvectors as a function of the normalized wavenumber in the form of leading terms in its series expansion and find that this error exceeds the dispersion error, except for the lowest degree where the eigenvector error is zero. We also present results of numerical experiments that confirm this analysis. Chapter 3 concluded that defect correction can improve the convergence property of finite-elements in the first-order system of acoustic equations in 1D; the inexpensive linear elements showed the same performance as a fourth-order scheme. However, for realistic problems we need to ensure that the same improvement holds in higher dimensions. Based on the results of the earlier chapter, we conjecture that defect correction should work for 2D problems. In the first half of chapter 4, we analyze the 2-D case. Theoretical results imply that the lowest-degree polynomial provides fourth-order accuracy with defect correction, if the grid of squares or triangles is highly regular and material properties constant. But numerical results converge more slowly than theoretical predictions. Further investigation demonstrates that this is due to the activation of error-inducing wavenumbers in the delta-source representation. In the second half of the chapter, we provide a solution to this problem in the form of a tapered-sinc source function. In chapter 5, we consider isotropic elastic wave propagation with continuous masslumped finite elements on tetrahedra with explicit time stepping. These elements require higher order polynomials in their interior to preserve accuracy after mass lumping and are recently discovered up to degree 4. Global assembly of the symmetric stiffness matrix is a natural approach but requires large memory. Local assembly on the fly, in the form of matrix-vector products per element at each time step, has a much smaller memory footprint. With dedicated expressions for local assembly, our code ran about 1.3 times faster for degree 2 and 1.9 times for degree 3 on a simple homogeneous test problem, using 24 cores. This is similar to the acoustic case. For a more realistic problem, the gain in efficiency was a factor 2.5 for degree 2 and 3 for degree 3. For the lowest degree, the linear element, the expressions for both the global and local assembly can be further simplified. In that case, global assembly is more efficient than local assembly. Among the three degrees, the element of degree 3 is the most efficient in terms of accuracy at a given cost. In chapter 6, we consider cubic Hermite elements as interpolants in place of Legendre polynomials. By nature of theirC1 continuity, they might offer a solution to the problems of ‘spurious’ wavenumbers seen in earlier chapters with conventional interpolation schemes. Results show acceptable convergence properties on homogeneous media, but the representation needs to be altered to suit discontinuities in density, which makes
interesting future work.