High-order space-time finite element schemes for acoustic and viscodynamic wave equations with temporal decoupling

Copyright @ 2014 The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.We revisit a method originally introduced by Werder et al. (in Comput. Methods Appl. Mech. Engrg., 190:6685–6708, 2001) for temporally discontinuous Galerkin FEMs applied to a parabolic partial differential equation. In that approach, block systems arise because of the coupling of the spatial systems through inner products of the temporal basis functions. If the spatial finite element space is of dimension D and polynomials of degree r are used in time, the block system has dimension (r + 1)D and is usually regarded as being too large when r > 1. Werder et al. found that the space-time coupling matrices are diagonalizable over inline image for r ⩽100, and this means that the time-coupled computations within a time step can actually be decoupled. By using either continuous Galerkin or spectral element methods in space, we apply this DG-in-time methodology, for the first time, to second-order wave equations including elastodynamics with and without Kelvin–Voigt and Maxwell–Zener viscoelasticity. An example set of numerical results is given to demonstrate the favourable effect on error and computational work of the moderately high-order (up to degree 7) temporal and spatio-temporal approximations, and we also touch on an application of this method to an ambitious problem related to the diagnosis of coronary artery disease

Optimally blended spectral-finite element scheme for wave propagation and nonstandard reduced integration

We study the dispersion and dissipation of the numerical scheme obtained by taking a weighted averaging of the consistent (finite element) mass matrix and lumped (spectral element) mass matrix for the small wave number limit. We find and prove that for the optimum blending the resulting scheme (a) provides $2p+4$ order accuracy for $p$th order method (two orders more accurate compared with finite and spectral element schemes); (b) has an absolute accuracy which is $\mathcal{O}(p^{-3})$ and $\mathcal{O}(p^{-2})$ times better than that of the pure finite and spectral element schemes, respectively; (c) tends to exhibit phase lag. Moreover, we show that the optimally blended scheme can be efficiently implemented merely by replacing the usual Gaussian quadrature rule used to assemble the mass and stiffness matrices by novel nonstandard quadrature rules which are also derived

Numerical Methods and Algorithms for High Frequency Wave Scattering Problems in Homogeneous and Random Media

This dissertation consists of four integral parts with a unified objective of developing efficient numerical methods for high frequency time-harmonic wave equations defined on both homogeneous and random media. The first part investigates the generalized weak coercivity of the acoustic Helmholtz, elastic Helmholtz, and time-harmonic Maxwell wave operators. We prove that such a weak coercivity holds for these wave operators on a class of more general domains called generalized star-shape domains. As a by-product, solution estimates for the corresponding Helmholtz-type problems are obtained. The second part of the dissertation develops an absolutely stable (i.e. stable in all mesh regimes) interior penalty discontinuous Galerkin (IP-DG) method for the elastic Helmholtz equations. A special mesh-dependent sesquilinear form is proposed and is shown to be weakly coercive in all mesh regimes. We prove that the proposed IP-DG method converges with optimal rate with respect to the mesh size. Numerical experiments are carried out to demonstrate the theoretical results and compare this method to the standard finite element method. The third part of the dissertation develops a Monte Carlo interior penalty discontinuous Galerkin (MCIP-DG) method for the acoustic Helmholtz equation defined on weakly random media. We prove that the solution to the random Helmholtz problem has a multi-modes expansion (i.e., a power series in a medium- related small parameter). Using this multi-modes expansion an efficient and accurate numerical method for computing moments of the solution to the random Helmholtz problem is proposed. The proposed method is also shown to converge optimally. Numerical experiments are carried out to compare the new multi-modes MCIP-DG method to a classical Monte Carlo method. The last part of the dissertation develops a theoretical framework for Schwarz pre- conditioning methods for general nonsymmetric and indefinite variational problems which are discretized by Galerkin-type discretization methods. Such a framework has been missing in the literature and is of great theoretical and practical importance for solving convection-diffusion equations and Helmholtz-type wave equations. Condition number estimates for the additive and hybrid Schwarz preconditioners are established under some structure assumptions. Numerical experiments are carried out to test the new framework

Seismic studies of 3-D elastic and anelastic structure of crust and upper mantle in Western China

Western China is a region with significant geological heterogeneity, especially because of the Tibetan Plateau, the largest and highest plateau in the world, which was created by the Cenozoic Indian-Eurasian continental collision. It is seismically very active even in some populous regions. Therefore, it is important to study the crustal and upper mantle structure beneath western China to understand the tectonics of the region and provide constraint for earthquake hazards. There are two important properties of Earth’s media: elastic and anelastic. Elastic properties, mainly seismic velocities, are studied by travel times of seismic waves. Anelastic properties are generally studied by attenuation of seismic waves retrieved from seismic amplitudes. We use different seismic datasets, including surface-wave dispersion from ambient noise and earthquakes, teleseismic receiver functions, and body-wave travel times, as well as joint inversions of these datasets utilizing neighborhood searching algorithm. We propose a generalized H-κ method with harmonic correction on receiver functions to provide a better estimation of the crustal thickness (H) and Vp/Vs ratio (κ). Our improved joint inversion results of S-wave velocities and crustal Vp/Vs ratios reveal mid-crustal low velocity, high Vp/Vs zones (possibly due to the presence of partial melt) in northern and southern Tibet, as well as in SE Tibet (in mid-lower crust), but not in central Tibet, which may suggest different characteristics of Lhasa Block in central Tibet. Uppermost mantle P and S velocities suggest that the subducted Indian mantle lithosphere is torn into at least 4 pieces that subduct at different angles and have different northern limits. This observation, when compared to seismicity and strain rates (from GPS), suggests coupled lithospheric deformation in southern Tibet. Notably, the lateral extent of potential megathrust earthquakes may be limited by the segment boundaries. The anelastic structure of the Earth, in particular seismic attenuation, is more useful in characterizing the temperature and fluid content as well as predicting earthquake ground motion. We have done simulations on the major factors that affect seismic amplitudes, such as scattering and focusing/defocusing. Strong amplification is observed within the major sedimentary basins, which sustains along the paths that pass through the thickest sediments and indicates that the source-to-basin direction of seismic waves affects their amplitudes. Internal scattering can generate strong coda which may interfere with the surface wave and make its amplitude difficult to measure. These factors need to be addressed with caution before the extraction of intrinsic attenuation structure in western China. Besides development of different seismic methods and models, this study makes considerable progress in answering two questions: Where earthquakes (especially large ones) are more likely to occur in western China, and what effects these earthquakes will have. Chapters 1-4 focus on the first issue by providing new seismic observations that can be compared with other geological constraints. This part includes the development of three seismic methods, namely generalized H-κ stacking method with harmonic corrections (Chapter 1), joint inversion of surface-wave dispersions and receiver functions with P-velocity model (Chapter 2), and joint inversion of dispersions, receiver functions, and Pn station delay time (Chapter 3), as well as comparison between seismic and geological observations and their implications (Chapter 4). Chapter 5 focuses on the second issue by numerically simulating the amplitudes of seismic waves, aiming at retrieving the attenuation structure in western China and predicting the amplitudes of future earthquakes