This paper deals with efficient numerical representation and manipulation of differential and integral operators as symbols in phase-space, i.e., functions of space @math and frequency @math . The symbol smoothness conditions obeyed by many operators in connection to smooth linear partial differential equations allow to write fast-converging, non-asymptotic expansions in adequate systems of rational Chebyshev functions or hierarchical splines. The classical results of closedness of such symbol classes under multiplication, inversion and taking the square root translate into practical iterative algorithms for realizing these operations directly in the proposed expansions. Because symbol-based numerical methods handle operators and not functions, their complexity depends on the desired resolution @math very weakly, typically only through @math factors. We present three applications to computational problems related to wave propagation: 1) preconditioning the Helmholtz equation, 2) decomposing wavefields into one-way components and 3) depth-stepping in reflection seismology.
The idea of writing pseudodifferential symbols in separated form to formulate various one-way approximations to the variable-coefficients Helmholtz equation has long been a tradition in seismic imaging. This almost invariably involves a high-frequency approximation of some kind. Some influential work includes the phase screen method by Fisk and McCartor @cite_34 , and the generalized screen expansion of Le Rousseau and de Hoop @cite_8 . This last reference discusses fast application of pseudodifferential operators in separated form using the FFT, and it is likely not the only reference to make this simple observation. A modern treatment of leading-order pseudodifferential approximations to one-way wave equations is in @cite_19 .
@cite_19 @cite_34 @cite_8
2113502771 2143055960 2083690159 2171798031
SUMMARY In this paper, we address the solution of three-dimensional heterogeneous Helmholtz problems discretized with second-order finite difference methods with application to acoustic waveform inversion in geophysics. In this setting, the numerical simulation of wave propagation phenomena requires the approximate solution of possibly very large indefinite linear systems of equations. For that purpose, we propose and analyse an iterative two-grid method acting on the Helmholtz operator where the coarse grid problem is solved inaccurately. A cycle of a multigrid method applied to a complex shifted Laplacian operator is used as a preconditioner for the approximate solution of this coarse problem. A single cycle of the new method is then used as a variable preconditioner of a flexible Krylov subspace method. We analyse the properties of the resulting preconditioned operator by Fourier analysis. Numerical results demonstrate the effectiveness of the algorithm on three-dimensional applications. The proposed numerical method allows us to solve three-dimensional wave propagation problems even at high frequencies on a reasonable number of cores of a distributed memory computer. Copyright © 2012 John Wiley & Sons, Ltd. In this work, we show how to incorporate attenuation into the optoacoustic tomography (OAT) imaging equation and develop a strategy for compensating for this attenuation during image reconstruction. In OAT, one exposes a sample to pulses of electromagnetic radiation that cause small amounts of heating in the specimen. The heating engenders thermal expansion which, in turn, gives rise to acoustic waves. The resulting acoustic pressure signal is generally measured by transducers arrayed around the object, and these data may be used to reconstruct images of the original electromagnetic absorption. Frequency-dependent absorption of the acoustic waves can lead to blurring and distortion in reconstructed images. We show that in the temporal frequency domain, the optoacoustic wave equation incorporating attenuation is equivalent to the inhomogeneous Helmholtz equation with a complex wave number. While some work has been done in other fields on directly solving Helmholtz equations with complex wave numbers, these are generally computationally intensive numerical approaches. We pursue a different approach, deriving an integral equation that relates the temporal optoacoustic signal at a given transducer location in the presence of attenuation to the ideal signal that would have been obtained in the absence of attenuation. This equation is readily discretized and the resulting linear system of equations involves a matrix that need only be inverted once, at which point the inverse can be used to correct all of the measured time signals prior to reconstruction by conventional methods. We present a novel boundary integral formulation of the Helmholtz transmission problem for bounded composite scatterers (that is, piecewise constant material parameters in "subdomains") that directly lends itself to operator preconditioning via Calderon projectors. The method relies on local traces on subdomains and weak enforcement of transmission conditions. The variational formulation is set in Cartesian products of standard Dirichlet and special Neumann trace spaces for which restriction and extension by zero are well defined. In particular, the Neumann trace spaces over each subdomain boundary are built as piecewise @math -distributions over each associated interface. Through the use of interior Calderon projectors, the problem is cast in variational Galerkin form with an operator matrix whose diagonal is composed of block boundary integral operators associated with the subdomains. We show existence and uniqueness of solutions based on an extension of Lions' projection lemma for non-closed subspaces. We also investigate asymptotic quasi-optimality of conforming boundary element Galerkin discretization. Numerical experiments in 2-D confirm the efficacy of the method and a performance matching that of another widely used boundary element discretization. They also demonstrate its amenability to different types of preconditioning. In some applications, it is reasonable to assume that geodesics (rays) have a consistent orientation so that the Helmholtz equation may be viewed as an evolution equation in one of the spatial directions. With such applications in mind, we propose a new Eulerian computational geometrical-optics method, dubbed the fast Huygens sweeping method, for computing Green functions of Helmholtz equations in inhomogeneous media in the high-frequency regime and in the presence of caustics. The first novelty of the new method is that the Huygens–Kirchhoff secondary source principle is used to integrate many locally valid asymptotic solutions to yield a globally valid asymptotic solution so that caustics associated with the usual geometrical-optics ansatz can be treated automatically. The second novelty is that a butterfly algorithm is adapted to carry out the matrix–vector products induced by the Huygens–Kirchhoff integration in O(NlogN) operations, where N is the total number of mesh points, and the proportionality constant depends on the desired accuracy and is independent of the frequency parameter. To reduce the storage of the resulting traveltime and amplitude tables, we compress each table into a linear combination of tensor-product based multivariate Chebyshev polynomials so that the information of each table is encoded into a small number of Chebyshev coefficients. The new method enjoys the following desired features: (1) it precomputes a set of local traveltime and amplitude tables; (2) it automatically takes care of caustics; (3) it constructs Green functions of the Helmholtz equation for arbitrary frequencies and for many point sources; (4) for a specified number of points per wavelength it constructs each Green function in nearly optimal complexity in terms of the total number of mesh points, where the prefactor of the complexity only depends on the specified accuracy and is independent of the frequency parameter. Both two-dimensional (2-D) and three-dimensional (3-D) numerical experiments are presented to demonstrate the performance and accuracy of the new method.