Source-linked AI summary
Sweeping Preconditioner for the Helmholtz Equation: Moving Perfectly Matched Layers
Björn Engquist, Lexing Ying
TL;DR
The paper targets the large, highly indefinite systems arising from high-frequency Helmholtz discretizations. It introduces a moving-PML sweeping preconditioner based on approximate layerwise LDLt factorization and lower-dimensional subproblem solves. The resulting method has linear application cost in 2D, near-linear application cost in 3D, and nearly frequency- and size-independent iteration counts in the reported tests.
Problem
High-frequency Helmholtz discretizations produce extremely large, highly indefinite systems with oscillatory Green’s functions that challenge standard direct and iterative solvers.
Method
The method approximates a layerwise LDLt factorization by representing Schur complements with interior moving PMLs, solved as quasi-1D banded-LU problems in 2D and quasi-2D multifrontal problems in 3D.
Results
The preconditioner has application costs O(N) in 2D and O(N log N) in 3D, while iteration counts are essentially independent of N or frequency; 3D tests required at most 13 GMRES iterations.
Takeaways & Limitations
Moving PMLs reduce the effective dimension of the subproblems, enabling fast direct solves that yield an effective sweeping preconditioner in both 2D and 3D.
Takeaways & Limitations
The study uses PML boundary conditions and basic central-difference stencils, whose dispersion accuracy is relatively poor; other boundary conditions and discretizations remain to be investigated.
Abstract
from arXiv · showhide
This paper introduces a new sweeping preconditioner for the iterative solution of the variable coefficient Helmholtz equation in two and three dimensions. The algorithms follow the general structure of constructing an approximate $LDL^t$ factorization by eliminating the unknowns layer by layer starting from an absorbing layer or boundary condition. The central idea of this paper is to approximate the Schur complement matrices of the factorization using moving perfectly matched layers (PMLs) introduced in the interior of the domain. Applying each Schur complement matrix is equivalent to solving a quasi-1D problem with a banded LU factorization in the 2D case and to solving a quasi-2D problem with a multifrontal method in the 3D case. The resulting preconditioner has linear application cost and the preconditioned iterative solver converges in a number of iterations that is essentially indefinite of the number of unknowns or the frequency. Numerical results are presented in both two and three dimensions to demonstrate the efficiency of this new preconditioner.
1 Introduction
The paper addresses the difficult, large-scale solution of high-frequency Helmholtz systems, where discretization produces O(ω^d) unknowns and oscillatory indefiniteness undermines standard solvers. It introduces moving-PML sweeping approximations that reduce subproblems by one dimension while retaining effective iterative performance.
- 1 Introduction: High-frequency discretization yields n proportional to ω per dimension, so the system has N = O(ω^d) unknowns and is extremely large.The Helmholtz system is also highly indefinite and has a very oscillatory Green’s function, making many multiscale-based solvers inefficient.
- 1 Introduction: The method constructs an approximate LDLt factorization by eliminating layers from an absorbing boundary and approximating intermediate Schur complements with moving PMLs.The moving PMLs are introduced in the domain interior and represent the Schur complement matrices through local Helmholtz solves.
- 1 Introduction: In 2D, moving-PML subproblems are quasi-1D systems solved by banded LU, while in 3D they are quasi-2D systems solved with multifrontal methods.This dimension reduction is intended to exploit efficient lower-dimensional direct solvers while producing an approximate inverse for iterative use.
- 1 Introduction: The preconditioner has application cost O(N) in 2D and O(N log N) in 3D, with construction costs O(N) and O(N4/3), respectively.These costs follow from the stated 2D and 3D moving-PML implementations.
2 Preconditioner in 2D
The 2D sweeping preconditioner approximates a layer-by-layer block LDL^t factorization by replacing dense Schur-complement applications with local moving-PML solves. These narrow subproblems can be applied rapidly, yielding an approximate inverse with nearly linear iterative-solver cost.
- Discretization: The 2D grid uses n points per dimension with N = n2 total points, where n is proportional to the wave number because each wavelength requires a constant number of samples.The discretization uses a Cartesian grid and a 5-point central-difference stencil.
- Sweeping factorization: The method begins with a block LDL^t factorization that eliminates grid layers from the absorbing boundary and represents each Schur complement through a half-space Green’s function.The discrete Helmholtz system is ordered row by row, and the factorization proceeds layer by layer from x2 = 0.
- Cost and approximation: The original dense sweeping factorization costs O(n4) = O(N2) to construct and O(n3) = O(N3/2) to apply, making it less useful than multifrontal methods by itself.These costs motivate replacing the dense Schur complements with moving-PML approximations.
- Moving PML: The moving-PML approximation pushes an artificial PML from x2 = 0 to immediately next to the active layer x2 = mh, restricting the auxiliary problem to a narrow subgrid.The moving PML is introduced only to approximate the Schur complement efficiently; it is not part of the original problem.
- Moving PML: Applying the approximate Schur complement requires solving a discrete Helmholtz system on a subgrid containing only b PML layers, which becomes a quasi-1D problem after reordering.The reordered matrix can be constructed efficiently, so applying the approximation is rapid.
- Cost and approximation: With moving PMLs, construction costs O(b3n2) = O(b3N), each preconditioner application costs O(n2) = O(N), and the iterative solver has nearly linear complexity when iterations grow at most logarithmically with N.The preconditioner is applied within GMRES, with total iterative cost O(NNI).
3 Numerical Results in 2D
The 2D experiments test the sweeping preconditioner across varied velocity fields, frequencies, discretization densities, and scattering settings. The results show nearly linear solution scaling with problem size and iteration counts that remain stable or weakly dependent on problem parameters.
- Varying frequency: When ω and n double, N quadruples, setup cost grows roughly fourfold, and total solution time grows by a factor of 4 or 5.The setup behavior is consistent with O(N) complexity, while the iteration count is essentially independent of n.
- Varying frequency: GMRES converges in about 20 iterations across the tested velocity fields with the sweeping preconditioner.The second and third velocity-field results exhibit quantitative behavior similar to the first.
- Varying discretization density: With ω/(2π) fixed at 32 and q varied from 8 to 64, iteration counts scale at most logarithmically and solution time scales roughly linearly with unknowns.These tests are summarized for all three velocity fields in Tables 4–6.
- Scattering problem: The 2D sweeping preconditioner is also applied to a sound-soft disk scattering problem by sweeping radially from the outer PML boundary toward the scatterer.The example uses c(r, θ) = 1, r0 = 0.15, r1 = 0.5, and q = 8 points per wavelength.
4 Preconditioner in 3D
The 3D preconditioner discretizes the Helmholtz equation on layers and approximates each Schur complement with a moving interior PML. These quasi-2D subproblems are solved by multifrontal methods within an approximate sweeping factorization.
- Discretization: The 3D domain is discretized on a Cartesian grid with n points per dimension, giving N = n^3 unknowns and a 7-point finite-difference operator.The grid spacing is h = 1/(n + 1), and n is proportional to the angular frequency.
- Sweeping factorization: The block system is tridiagonal by layers, with n^2 × n^2 diagonal blocks and diagonal off-diagonal blocks.The unknowns and right-hand side are grouped by the third coordinate, producing the layer vectors used by the sweep.
- Moving PML approximation: Each Schur complement maps a force on one layer to the solution restricted to that layer, and is approximated by moving a PML next to the current layer.The auxiliary domain contains b layers, yielding a localized approximation to the layer Schur complement.
- Moving PML approximation: The resulting operator is a quasi-2D Helmholtz problem solved efficiently with a modified multifrontal method.The method hierarchically partitions the (x1, x2) plane while retaining the layer direction in each group.
- Approximate inversion and preconditioner: The modified operator Aα is used for stability, and GMRES applies its approximate inverse as a preconditioner with total cost O(N_I N log N).The iteration count N_I is reported as essentially independent of N, giving almost linear overall complexity.
5 Numerical Results in 3D
Three-dimensional tests use varied velocity fields, external forces, frequencies, and discretization densities. The preconditioned GMRES iteration count stays nearly constant, while solution cost scales roughly linearly with the number of unknowns.
- Test problems: The experiments use a converging lens, a vertical waveguide, and a random velocity field, each tested with two external forces.The forces are a Gaussian point source and a Gaussian wave packet, producing waves that bend, form caustics, or propagate as a beam.
- Varying frequency: The frequency tests use q = 8 points per wavelength, n = 39, 79, 159, PML width b = 6, and three layers processed per iteration.The two sweeping fronts start from x3 = 0 and x3 = 1.
- Varying frequency: At most 12 preconditioned GMRES iterations are required for the first velocity field across the tested frequencies.For ω/(2π) = 20, the problem has more than four million unknowns and takes less than 600 seconds to solve.
- Varying frequency: Across all three velocity fields, the frequency tests require at most 13 GMRES iterations with the new sweeping preconditioner.The reported results cover the first, second, and third velocity fields.
- Varying discretization density: When q varies from 8 to 32 at fixed ω/(2π) = 5, iteration counts remain roughly constant and solution time scales roughly linearly with unknowns.This behavior is reported for all three velocity fields.
- Comparison with previous work: For a 20-wavelength problem, the current setup time is in the hundreds of seconds, compared with tens of thousands of seconds for the previous 3D implementation.The paper attributes the reduction mainly to a more efficient multifrontal implementation.
6 Conclusion and Future Work
The paper concludes that moving PMLs provide a flexible sweeping preconditioner for two- and three-dimensional Helmholtz problems. It identifies more accurate discretizations and scalable parallel implementations as future directions.
- Conclusion: The proposed method combines layer-by-layer approximate block LDLt elimination with moving interior PMLs to approximate Schur complements in two and three dimensions.The 2D subproblems use optimally ordered banded LU, while the 3D subproblems use multifrontal methods.
- Conclusion: The preconditioner has linear application cost, and GMRES iterations are essentially independent of the number of unknowns or frequency.These properties are the principal reported algorithmic outcomes of the construction.
- Future work: The current study uses PML boundary conditions and basic central-difference stencils, leaving more accurate stencils and finite-element-type discretizations for future investigation.The paper specifically mentions finite element, spectral element, and discontinuous Galerkin discretizations.
- Future work: Large-scale 3D use requires parallel processing, although multifrontal calculations within each sweep can be well parallelized.The paper also identifies coarse-scale domain decomposition as an alternative parallelization strategy.
- Flexibility: The approach is readily applicable to non-uniform and adaptive grids when those grids resolve the moving PMLs sufficiently accurately.Available multifrontal methods for such grids support this flexibility.