We present a solver for the 2D high-frequency Helmholtz equation in heterogeneous acoustic media, with online parallel complexity that scales optimally as @math , where @math is the number of volume unknowns, and @math is the number of processors, as long as @math grows at most like a small fractional power of @math . The solver decomposes the domain into layers, and uses transmission conditions in boundary integral form to explicitly define "polarized traces", i.e., up- and down-going waves sampled at interfaces. Local direct solvers are used in each layer to precompute traces of local Green's functions in an embarrassingly parallel way (the offline part), and incomplete Green's formulas are used to propagate interface data in a sweeping fashion, as a preconditioner inside a GMRES loop (the online part). Adaptive low-rank partitioning of the integral kernels is used to speed up their application to interface data. The method uses second-order finite differences. The complexity scalings are empirical but motivated by an analysis of ranks of off-diagonal blocks of oscillatory integrals. They continue to hold in the context of standard geophysical community models such as BP and Marmousi 2, where convergence occurs in 5 to 10 GMRES iterations.
Bhowmik and Stolk recently proposed new rules of optimized coarse grid corrections for a two level multigrid method @cite_65 ; Brandt and Livshits developed the wave-ray method @cite_68 , in which the oscillatory error components are eliminated by a ray-cycle, that exploits a geometric optics approximation of the Green's function; Elman and Ernst use a relaxed multigrid method as a preconditioner for an outer Krylov iteration in @cite_53 ; Haber and McLachlan proposed an alternative formulation for the Hemholtz equation @cite_45 , which reduces the problem to solve an eikonal equation and an advection-diffusion-reaction equation, which can be solved efficiently by a Krylov method using a multigrid method as a preconditioner; and Grote and Schenk @cite_50 proposed an algebraic multi-level preconditioner for the Helmholtz equation.
@cite_53 @cite_65 @cite_45 @cite_50 @cite_68
2151478420 2001099567 2025033760 175176216
An iterative solution method, in the form of a preconditioner for a Krylov subspace method, is presented for the Helmholtz equation. The preconditioner is based on a Helmholtz-type differential operator with a complex term. A multigrid iteration is used for approximately inverting the preconditioner. The choice of multigrid components for the corresponding preconditioning matrix with a complex diagonal is validated with Fourier analysis. Multigrid analysis results are verified by numerical experiments. High wavenumber Helmholtz problems in heterogeneous media are solved indicating the performance of the preconditioner. 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. An algebraic multilevel (ML) preconditioner is presented for the Helmholtz equation in heterogeneous media. It is based on a multilevel incomplete @math factorization and preserves the inherent (complex) symmetry of the Helmholtz equation. The ML preconditioner incorporates two key components for efficiency and numerical stability: symmetric maximum weight matchings and an inverse-based pivoting strategy. The former increases the block-diagonal dominance of the system, whereas the latter controls @math for numerical stability. When applied recursively, their combined effect yields an algebraic coarsening strategy, similar to algebraic multigrid methods, even for highly indefinite matrices. The ML preconditioner is combined with a Krylov subspace method and applied as a “black-box” solver to a series of challenging two- and three-dimensional test problems, mainly from geophysical seismic imaging. The numerical results demonstrate the robustness and efficiency of the ML preconditioner, even at higher frequency regimes. Using a finite element method to solve the Helmholtz equation leads to a sparse system of equations which in three dimensions is too large to solve directly. It is also non-Hermitian and highly indefinite and consequently difficult to solve iteratively. The approach taken in this paper is to precondition this linear system with a new preconditioner and then solve it iteratively using a Krylov subspace method. Numerical analysis shows the preconditioner to be effective on a simple 1D test problem, and results are presented showing considerable convergence acceleration for a number of different Krylov methods for more complex problems in 2D, as well as for the more general problem of harmonic disturbances to a non-stagnant steady flow. In this paper, we consider the numerical solution of the Helmholtz equation, arising from the study of the wave equation in the frequency domain. The approach proposed here differs from those recently considered in the literature, in that it is based on a decomposition that is exact when considered analytically, so the only degradation in computational performance is due to discretization and roundoff errors. In particular, we make use of a multiplicative decomposition of the solution of the Helmholtz equation into an analytical plane wave and a multiplier, which is the solution of a complex-valued advection-diffusion-reaction equation. The use of fast multigrid methods for the solution of this equation is investigated. Numerical results show that this is an efficient solution algorithm for a reasonable range of frequencies. In this paper we survey the development of fast iterative solvers aimed at solving 2D 3D Helmholtz problems. In the first half of the paper, a survey on some recently developed methods is given. The second half of the paper focuses on the development of the shifted Laplacian preconditioner used to accelerate the convergence of Krylov subspace methods applied to the Helmholtz equation. Numerical examples are given for some difficult problems, which had not been solved iteratively before. List of Algorithms Preface 1. Introduction. Brief Overview of the State of the Art Notation Review of Relevant Linear Algebra Part I. Krylov Subspace Approximations. 2. Some Iteration Methods. Simple Iteration Orthomin(1) and Steepest Descent Orthomin(2) and CG Orthodir, MINRES, and GMRES Derivation of MINRES and CG from the Lanczos Algorithm 3. Error Bounds for CG, MINRES, and GMRES. Hermitian Problems-CG and MINRES Non-Hermitian Problems-GMRES 4. Effects of Finite Precision Arithmetic. Some Numerical Examples The Lanczos Algorithm A Hypothetical MINRES CG Implementation A Matrix Completion Problem Orthogonal Polynomials 5. BiCG and Related Methods. The Two-Sided Lanczos Algorithm The Biconjugate Gradient Algorithm The Quasi-Minimal Residual Algorithm Relation Between BiCG and QMR The Conjugate Gradient Squared Algorithm The BiCGSTAB Algorithm Which Method Should I Use? 6. Is There A Short Recurrence for a Near-Optimal Approximation? The Faber and Manteuffel Result Implications 7. Miscellaneous Issues. Symmetrizing the Problem Error Estimation and Stopping Criteria Attainable Accuracy Multiple Right-Hand Sides and Block Methods Computer Implementation Part II. Preconditioners. 8. Overview and Preconditioned Algorithms. 9. Two Example Problems. The Diffusion Equation The Transport Equation 10. Comparison of Preconditioners. Jacobi, Gauss--Seidel, SOR The Perron--Frobenius Theorem Comparison of Regular Splittings Regular Splittings Used with the CG Algorithm Optimal Diagonal and Block Diagonal Preconditioners 11. Incomplete Decompositions. Incomplete Cholesky Decomposition Modified Incomplete Cholesky Decomposition 12. Multigrid and Domain Decomposition Methods. Multigrid Methods Basic Ideas of Domain Decomposition Methods. A preconditioned iterative solver for the 3D frequency-domain wave equation applied to seismic problems is evaluated. The preconditioner corresponds to an approximate inverse of a heavily damped wave equation deduced from the (undamped) wave equation. The approximate inverse is computed with one multigrid cycle. Numerical results show that the method is robust and that the number of iterations increases roughly linearly with frequency when the grid spacing is adapted to keep a constant number of discretization points per wavelength. To evaluate the relevance of this iterative solver, the number of floating-point operations required for two imaging problems are roughly evaluated. This rough estimate indicates that the time-domain migration approach is more than one order of magnitude faster. The full-wave-form tomography, based on a least-squares formulation and a scale separation approach, has the same complexity in both domains. We present an incomplete LU preconditioner for solving discretized Helmholtz problems. The preconditioner is based on an analytic factorization of the Helmholtz operator. This allows us to take the physical properties of the acoustics problem modeled by the Helmholtz equation into account in the preconditioner. We show how the parameters in the preconditioner can be chosen in order to make it effective. Numerical experiments show that the new preconditioner leads to convergent iterative methods even for large wave numbers, and it outperforms classical ILU preconditioners by a large margin. Standard multigrid algorithms have proven ineffective for the solution of discretizations of Helmholtz equations. In this work we modify the standard algorithm by adding GMRES iterations at coarse levels and as an outer iteration. We demonstrate the algorithm's effectiveness through theoretical analysis of a model problem and experimental results. In particular, we show that the combined use of GMRES as a smoother and outer iteration produces an algorithm whose performance depends relatively mildly on wave number and is robust for normalized wave numbers as large as 200. For fixed wave numbers, it displays grid-independent convergence rates and has costs proportional to the number of unknowns. In 1983, a preconditioner was proposed [J. Comput. Phys. 49 (1983) 443] based on the Laplace operator for solving the discrete Helmholtz equation efficiently with CGNR. The preconditioner is especially effective for low wavenumber cases where the linear system is slightly indefinite. Laird [Preconditioned iterative solution of the 2D Helmholtz equation, First Year's Report, St. Hugh's College, Oxford, 2001] proposed a preconditioner where an extra term is added to the Laplace operator. This term is similar to the zeroth order term in the Helmholtz equation but with reversed sign. In this paper, both approaches are further generalized to a new class of preconditioners, the so-called "shifted Laplace" preconditioners of the form Δφ-αk2φ with α ∈ C. Numerical experiments for various wavenumbers indicate the effectiveness of the preconditioner. The preconditioner is evaluated in combination with GMRES, Bi-CGSTAB, and CGNR. Linear systems which originate from the simulation of wave propagation phenomena can be very difficult to solve by iterative methods. These systems are typically complex valued and they tend to be highly indefinite, which renders the standard ILU-based preconditioners ineffective. This paper presents a study of ways to enhance standard preconditioners by altering the diagonal by imaginary shifts. Prior work indicates that modifying the diagonal entries during the incomplete factorization process, by adding to it purely imaginary values can improve the quality of the preconditioner in a substantial way. Here we propose simple algebraic heuristics to perform the shifting and test these techniques with the ARMS and ILUT preconditioners. Comparisons are made with applications stemming from the diffraction of an acoustic wave incident on a bounded obstacle (governed by the Helmholtz Wave Equation). We investigate the parallel performance of an iterative solver for 3D heterogeneous Helmholtz problems related to applications in seismic wave propagation. For large 3D problems, the computation is no longer feasible on a single processor, and the memory requirements increase rapidly. Therefore, parallelization of the solver is needed. We employ a complex shifted-Laplace preconditioner combined with the Bi-CGSTAB iterative method and use a multigrid method to approximate the inverse of the resulting preconditioning operator. A 3D multigrid method with 2D semi-coarsening is employed. We show numerical results for large problems arising in geophysical applications. The Helmholtz equation governing wave propagation and scattering phenomena is difficult to solve numerically. Its discretization with piecewise linear finite elements results in typically large linear systems of equations. The inherently parallel domain decomposition methods constitute hence a promising class of preconditioners. An essential element of these methods is a good coarse space. Here, the Helmholtz equation presents a particular challenge, as even slight deviations from the optimal choice can be devastating. In this paper, we present a coarse space that is based on local eigenproblems involving the Dirichlet-to-Neumann operator. Our construction is completely automatic, ensuring good convergence rates without the need for parameter tuning. Moreover, it naturally respects local variations in the wave number and is hence suited also for heterogeneous Helmholtz problems. The resulting method is parallel by design and its efficiency is demonstrated on 2D homogeneous and heterogeneous numerical examples. We propose and study a domain decomposition method for solving the truncated perfectly matched layer (PML) approximation in bounded domain of Helmholtz scattering problems. The method is based on the decomposition of the domain into nonoverlapping layers and the idea of source transfer which transfers the sources equivalently layer by layer so that the solution in the final layer can be solved using a PML method defined locally outside the last two layers. The convergence of the method is proved for the case of constant wave number based on the analysis of the fundamental solution of the PML equation. The method can be used as an efficient preconditioner in the preconditioned GMRES method for solving discrete Helmholtz equations with constant and heterogeneous wave numbers. Numerical examples are included. The paper introduces the sweeping preconditioner, which is highly efficient for iterative solutions of the variable-coefficient Helmholtz equation including very-high-frequency problems. The first central idea of this novel approach is to construct an approximate factorization of the discretized Helmholtz equation by sweeping the domain layer by layer, starting from an absorbing layer or boundary condition. Given this specific order of factorization, the second central idea is to represent the intermediate matrices in the hierarchical matrix framework. In two dimensions, both the construction and the application of the preconditioners are of linear complexity. The generalized minimal residual method (GMRES) solver with the resulting preconditioner converges in an amazingly small number of iterations, which is essentially independent of the number of unknowns. This approach is also extended to the three-dimensional case with some success. Numerical results are provided in both two and three dimensions to demonstrate the efficiency of this new approach. These notes were written for an introductory course on the application of multigrid methods to elliptic and hyperbolic partial differential equations for engineers, physicists and applied mathematicians. The use of more advanced mathematical tools, such as functional analysis, is avoided. The course is intended to be accessible to a wide audience of users of computational methods. We restrict ourselves to finite volume and finite difference discretization. The basic principles are given. Smoothing methods and Fourier smoothing analysis are reviewed. The fundamental multigrid algorithm is studied. The smoothing and coarse grid approximation properties are discussed. Multigrid schedules and structured programming of multigrid algorithms are treated. Robustness and efficiency are considered. We describe new algorithms of the locally optimal block preconditioned conjugate gradient (LOBPCG) method for symmetric eigenvalue problems, based on a local optimization of a three-term recurrence, and suggest several other new methods. To be able to compare numerically different methods in the class, with different preconditioners, we propose a common system of model tests, using random preconditioners and initial guesses. As the "ideal" control algorithm, we advocate the standard preconditioned conjugate gradient method for finding an eigenvector as an element of the null-space of the corresponding homogeneous system of linear equations under the assumption that the eigenvalue is known. We recommend that every new preconditioned eigensolver be compared with this "ideal" algorithm on our model test problems in terms of the speed of convergence, costs of every iteration, and memory requirements. We provide such comparison for our LOBPCG method. Numerical results establish that our algorithm is practically as efficient as the ideal'' algorithm when the same preconditioner is used in both methods. We also show numerically that the LOBPCG method provides approximations to first eigenpairs of about the same quality as those by the much more expensive global optimization method on the same generalized block Krylov subspace. We propose a new version of block Davidson's method as a generalization of the LOBPCG method. Finally, direct numerical comparisons with the Jacobi--Davidson method show that our method is more robust and converges almost two times faster. Abstract Efficient algorithms for the solution of partial differential equations on parallel computers are often based on domain decomposition methods. Schwarz preconditioners combined with standard Krylov space solvers are widely used in this context, and such a combination is shown here to perform very well in the case of the Wilson–Dirac equation in lattice QCD. In particular, with respect to even-odd preconditioned solvers, the communication overhead is significantly reduced, which allows the computational work to be distributed over a large number of processors with only small parallelization losses. We consider optimal-scaling multigrid solvers for the linear systems that arise from the discretization of problems with evolutionary behavior. Typically, solution algorithms for evolu- tion equations are based on a time-marching approach, solving sequentially for one time step after the other. Parallelism in these traditional time-integration techniques is limited to spatial parallelism. However, current trends in computer architectures are leading toward systems with more, but not faster, processors. Therefore, faster compute speeds must come from greater parallelism. One ap- proach to achieving parallelism in time is with multigrid, but extending classical multigrid methods for elliptic operators to this setting is not straightforward. In this paper, we present a nonintrusive, optimal-scaling time-parallel method based on multigrid reduction (MGR). We demonstrate optimal- ity of our multigrid-reduction-in-time algorithm (MGRIT) for solving diffusion equations in two and three space dimensions in numerical experiments. Furthermore, through both parallel performance models and actual parallel numerical results, we show that we can achieve significant speedup in comparison to sequential time marching on modern architectures.