2D Helmholtz Equation Research Articles

The method of polarized traces for the 3D Helmholtz equation

We have developed a fast solver for the 3D Helmholtz equation, in heterogeneous, constant density, acoustic media, in the high-frequency regime. The solver is based on the method of polarized traces, a layered domain-decomposition method, where the subdomains are connected via transmission conditions prescribed by the discrete Green’s representation formula and artificial reflections are avoided by enforcing nonreflecting boundary conditions between layers. The method of polarized traces allows us to consider only unknowns at the layer interfaces, reducing the overall cost and memory footprint of the solver. We determine that polarizing the wavefields in this manner yields an efficient preconditioner for the reduced system, whose rate of convergence is independent of the problem frequency. The resulting preconditioned system is solved iteratively using generalized minimum residual, where we never assemble the reduced system or preconditioner; rather, we implement them via solving the Helmholtz equation locally within the subdomains. The method is parallelized using Message Passing Interface and coupled with a distributed linear algebra library and pipelining to obtain an empirical on-line runtime [Formula: see text], where [Formula: see text] is the total number of degrees of freedom, [Formula: see text] is the number of subdomains, and [Formula: see text] is the number of right-hand sides (RHS). This scaling is favorable for regimes in which the number of sources (distinct RHS) is large, for example, enabling large-scale implementations of frequency-domain full-waveform inversion.

The modified method of fundamental solutions for exterior problems of the Helmholtz equation; spurious eigenvalues and their removals

For the exterior Dirichlet problems (EDP) of the Helmholtz equation in 2D, when the Sommerfeld radiation condition is satisfied, there always exists the unique solution. Denote the unbounded domain S∞ outside of a bounded simply-connected domain, where the interior boundary Γin is non-circular. The standard MFS is first studied by using the complex Hankel function H0(1)(kr) as the fundamental solutions (FS). However, some solutions can not be obtained correctly (e.g., the spurious eigenvalues called). This paper is devoted to explore the spurious eigenvalues and their removals. For the non-circular interior boundary Γin, the error analysis is briefly made for the standard MFS, and some guidance is provided to bypass the spurious eigenvalues. Denote two nodes: P=(ρ,θ) and Q=(R,ϕ), and r=|PQ‾|=R2+ρ2−2Rρcos⁡(θ−ϕ). To completely eliminate all spurious eigenvalues, we solicit the new combined FS: (∂∂R±ik)H0(1)(kr) with i=−1, and propose the new modified MFS. The new algorithms are simple, and the strict analysis may be made. The bounds of errors are derived, and the polynomial convergence rates can be achieved. The bounds of condition numbers are derived for circular Γ∘in only, to display the exponential growth via the number of the new FS used. Numerical experiments are carried to support the analysis made. This paper is the first time for the MFS to explore the spurious eigenvalues and their removals, accompanied with strict analysis of errors and stability.

A geometric full multigrid method and fourth-order compact scheme for the 3D Helmholtz equation on nonuniform grid discretization

A full multigrid method and fourth-order compact difference scheme is designed to solve the 3D Helmholtz equation on unequal mesh size. Three dimensional restriction and prolongation operators of the multigrid method on unequal grids could be constructed based on volume law. Two numerical experiments are implemented, and the results show the computational efficiency and accuracy of the full multigrid method. The study also illustrates that the full multigrid method with fourth-order compact difference scheme has great advantages in computation, which has been time consuming, and in iterative convergence efficiency.

Radial basis function-generated finite differences with Bessel weights for the 2D Helmholtz equation

In this paper we obtain approximated numerical solutions for the 2D Helmholtz equation using a Radial Basis Function-generated Finite Difference (RBF-FD) scheme, where weights are calculated by applying an oscillatory radial basis function given in terms of Bessel functions of the first kind. The problem of obtaining weights by local interpolation is ill-conditioned; we overcome this difficulty by means of regularization of the interpolation matrix by perturbing its diagonal. The condition number of this perturbed matrix is controlled according to a prescribed value of a regularization parameter. Different numerical tests are performed in order to study convergence and algorithmic complexity. As a result, we verify that dispersion and pollution effects are mitigated.

A Unified 2D/3D Large-Scale Software Environment for Nonlinear Inverse Problems

Large-scale parameter estimation problems are among some of the most computationally demanding problems in numerical analysis. An academic researcher’s domain-specific knowledge often precludes that of software design, which results in inversion frameworks that are technically correct but not scalable to realistically sized problems. On the other hand, the computational demands for realistic problems result in industrial codebases that are geared solely for high performance, rather than comprehensibility or flexibility. We propose a new software design for inverse problems constrained by partial differential equations that bridges the gap between these two seemingly disparate worlds. A hierarchical and modular design reduces the cognitive burden on the user while exploiting high-performance primitives at the lower levels. Our code has the added benefit of actually reflecting the underlying mathematics of the problem, which lowers the cognitive load on the user using it and reduces the initial startup period before a researcher can be fully productive. We also introduce a new preconditioner for the 3D Helmholtz equation that is suitable for fault-tolerant distributed systems. Numerical experiments on a variety of 2D and 3D test problems demonstrate the effectiveness of this approach on scaling algorithms from small- to large-scale problems with minimal code changes.

Robust Adaptive $hp$ Discontinuous Galerkin Finite Element Methods for the Helmholtz Equation

This paper presents an $hp$ a posteriori error analysis for the 2D Helmholtz equation that is robust in the polynomial degree $p$ and the wave number $k$. For the discretization, we consider a discontinuous Galerkin formulation that is unconditionally well posed. The a posteriori error analysis is based on the technique of equilibrated fluxes applied to a shifted Poisson problem, with the error due to the nonconformity of the discretization controlled by a potential reconstruction. We prove that the error estimator is both reliable and efficient, under the condition that the initial mesh size and polynomial degree are chosen such that the discontinuous Galerkin formulation converges, i.e., it is out of the regime of pollution. We confirm the efficiency of an $hp$-adaptive refinement strategy based on the presented robust a posteriori error estimator via several numerical examples.

A fast singular boundary method for 3D Helmholtz equation

This paper presents a fast singular boundary method (SBM) for three-dimensional (3D) Helmholtz equation. The SBM is a boundary-type meshless method which incorporates the advantages of the boundary element method (BEM) and the method of fundamental solutions (MFS). It is easy-to-program, and attractive to the problems with complex geometries. However, the SBM is usually limited to small-scale problems, because of the operation count of O(N3) with direct solvers or O(N2) with iterative solvers, as well as the memory requirement of O(N2). To overcome this drawback, this study makes the first attempt to employ the precorrected-FFT (PFFT) to accelerate the SBM matrix–vector multiplication at each iteration step of the GMRES for 3D Helmholtz equation. Consequently, the computational complexity can be reduced from O(N2) to O(NlogN) or O(N). Three numerical examples are successfully tested on a desktop computer. The results clearly demonstrate the accuracy and efficiency of the developed fast PFFT-SBM strategy.

New formulations for the hybrid meshless‐MoM method applied to 2D scattering problems

AbstractIn this work, new formulations for the hybrid meshless‐method of moments technique are presented. The proposed approaches are applied for solving 2D electromagnetic scattering problems. The solution is obtained by dividing the problem into 2 subproblems: one internal to the scatterer and another external to it. For the internal subproblem, the meshless method is used for solving the weak form of the 2D Helmholtz equation. For the external subproblem, the Method of Moments is used for solving different combinations of electric field integral equation and magnetic field integral equation. Both subproblems are coupled by imposing the continuity of the tangential fields over the scatterer surface. The technique is applied for solving the TMz electromagnetic scattering by an infinite dielectric cylinder. Numerical and analytical results are compared to each other showing the validity and accuracy of the proposed approaches.

A fourth-order orthogonal spline collocation solution to 1D-Helmholtz equation with discontinuity

In this paper, we use orthogonal spline collocation methods (OSCM) for the one dimensional Helmholtz equation with discontinuous coefficients. We discuss the existence uniqueness results and establish optimal error estimates. We use piecewise Hermite cubic basis functions to approximate the solution. Finally, we perform some numerical experiments and validate the theoretical results. Comparative to existing methods, we prove that the orthogonal spline collocation methods (OSCM) handles the discontinuous coefficients effectively and gives optimal order of convergence for $$\Vert y-y_h\Vert _{L^{\infty }}$$ -norm and superconvergent result for $$\Vert y'-y'_h\Vert _{L^{\infty }}$$ -norm at the grid points.

A framework of shape optimisation based on the isogeometric boundary element method toward designing thin-silicon photovoltaic devices

We propose a gradient-based optimisation framework to design the shape of the interfaces in periodic layered structures, which can model photovoltaic devices, optical gratings, etc., according to a desired electromagnetic property. To this end, we first develop an isogeometric boundary element method (IGBEM) to analyse singly periodic boundary value problems for 2D Helmholtz equation in layered media. By nature, the IGBEM is not only numerically accurate but also suitable for the shape optimisation in comparison with the conventional boundary- and domain-type solvers. Next, we derive the shape derivative (or shape sensitivity) of the objective function in terms of the magnetic field strength or energy absorption rate on the basis of the adjoint variable method. The shape derivative can be computed by solving the primal and adjoint problems with the IGBEM. Subsequently, we map our shape optimisation problems to nonlinear programming problems in order to exploit a general-purpose software for the latter problems. After verifying our optimisation framework rigorously by comparing with the exact solutions, we demonstrate a shape optimisation of the silicon-metal interface in a model of thin-silicon photovoltaic devices: it is crucial to determine the interface so that the energy of the incident light can be confined in the silicon layer as much as possible when the thickness of the layer is unconventionally small; one micrometre in our model. The optimal shape achieved 8.6 times higher absorption rate than the reference (flat) shape. This simulation shows that the proposed framework is capable of making a fundamental contribution to analysing and designing photovoltaic devices as well as photonic and plasmonic devices that consist of singly periodic layers.

Convexification of a 3-D coefficient inverse scattering problem

A version of the so-called “convexification” numerical method for a coefficient inverse scattering problem for the 3D Helmholtz equation is developed analytically and tested numerically. Backscattering data are used, which result from a single direction of the propagation of the incident plane wave on an interval of frequencies. The method converges globally. The idea is to construct a weighted Tikhonov-like functional. The key element of this functional is the presence of the so-called Carleman Weight Function (CWF). This is the function which is involved in the Carleman estimate for the Laplace operator. This functional is strictly convex on any appropriate ball in a Hilbert space for an appropriate choice of the parameters of the CWF. Thus, both the absence of local minima and convergence of minimizers to the exact solution are guaranteed. Numerical tests demonstrate a good performance of the resulting algorithm. Unlikeprevious the so-called tail functions globally convergent method, we neither do not impose the smallness assumption of the interval of wavenumbers, nor we do not iterate with respect to the so-called tail functions.

A Fast Fourth-Order Method for 3D Helmholtz Equation with Neumann Boundary

We present fast fourth-order finite difference scheme for 3D Helmholtz equation with Neumann boundary condition. We employ the discrete Fourier transform operator and divide the problem into some independent subproblems. By means of the Gaussian elimination in the vertical direction, the problem is reduced into a small system on the top layer of the domain. The procedure for solving the numerical solutions is accelerated by the sparsity of Fourier operator under the space complexity of O(M3). Furthermore, the method makes it possible to solve the 3D Helmholtz equation with large grid number. The accuracy and efficiency of the method are validated by two test examples which have exact solutions.

Nested Domain Decomposition with Polarized Traces for the 2D Helmholtz Equation

We present a solver for the two-dimensional high-frequency Helmholtz equation in heterogeneous, constant density, acoustic media, with online parallel complexity that scales empirically as ${\cal O}(\frac{N}{P})$, where $N$ is the number of volume unknowns, and $P$ is the number of processors, as long as $P = {\cal O}(N^{1/5})$. This sublinear scaling is achieved by domain decomposition, not distributed linear algebra, and improves on the $P = {\cal O}(N^{1/8})$ scaling reported earlier in [L. Zepeda-Nún͂ez and L. Demanet, J. Comput. Phys., 308 (2016), pp. 347--388]. The solver relies on a two-level nested domain decomposition: a layered partition on the outer level and a further decomposition of each layer in cells at the inner level. The Helmholtz equation is reduced to a surface integral equation (SIE) posed at the interfaces between layers, efficiently solved via a nested version of the polarized traces preconditioner [L. Zepeda-Nún͂ez and L. Demanet, J. Comput. Phys., 308 (2016), pp. 347--388]. The favorable complexity is achieved via an efficient application of the integral operators involved in the SIE.

A Fast High Order Algorithm for 3D Helmholtz Equation with Dirichlet Boundary

Helmholtz equation is widely applied in the scientific and engineering problem. For the solution of the three-dimensional Helmholtz equation, the computational efficiency of the algorithm is especially important. In this paper, in order to solve the contradiction between accuracy and efficiency, a fast high order finite difference method is proposed for solving the three-dimensional Helmholtz equation. First, a traditional fourth order method is constructed. Then, fast Fourier transformation are introduced to generate a block-tridiagonal structure which can easily divide the original problem into small and independent subsystems. For large 3D problems, the computation of traditional discrete Fourier transformation is less efficient, and the memory requirements increase rapidly. To fix this problem, the transformed coefficient matrix is constructed as a sparse structure. In light of the sparsity, the algorithm presented in this paper requires less memory space and computational cost. This sparse structure also leads to independent solving procedure of any plane in the domain. Therefore, parallel method can be used to solve the Helmholtz equation with large grid number. In the end, three numerical experiments are presented to verify the effectiveness of the fast fourth-order algorithm, and the acceleration effect to use the parallel method has been demonstrated.

Accelerated over relaxation iterative method using triangle element approximation for solving 2D Helmholtz Equations

Weighted iterative methods particularly Accelerated Over Relaxation (AOR) method are used to solve linear system generated from triangle finite element approximation equation in solving 2D Helmholtz equation. The development of the AOR iterative method were also presented. Numerical experiments have been carried out and the results obtained confirm the superiority of the proposed iterative method

High Order Compact Scheme and the Moving Mesh Method for Helmholtz Equation with Variable Wave Numbers

High order compact scheme has extensively been used for solving the Helmholtz equation on uniform mesh. Although uniform spacing is computationally advantageous for isotropic problems, the non-uniform mesh will be superior to the problems with the anisotropic properties or boundary layers. In this paper, we proposed a compact scheme on nonuniform mesh which has at least third order accuracy for the 2D Helmholtz equation with variable wave number. A novel adaptive mesh algorithm is also presented to generate the orthogonal nonuniform mesh. The 2D Helmholtz equation is solved by combining the HOC scheme on nonuniform mesh with the moving mesh method. The multigrid method is used to solve the discrete algebraic system. Numerical experiments are executed to verify the accuracy and efficiency of the presented method.

An iterative solver for the 3D Helmholtz equation

We develop a frequency–domain iterative solver for numerical simulation of acoustic waves in 3D heterogeneous media. It is based on the application of a unique preconditioner to the Helmholtz equation that ensures convergence for Krylov subspace iteration methods. Effective inversion of the preconditioner involves the Fast Fourier Transform (FFT) and numerical solution of a series of boundary value problems for ordinary differential equations. Matrix-by-vector multiplication for iterative inversion of the preconditioned matrix involves inversion of the preconditioner and pointwise multiplication of grid functions. Our solver has been verified by benchmarking against exact solutions and a time-domain solver.

Optimization of plane wave directions in plane wave discontinuous Galerkin methods for the Helmholtz equation

Recently, the use of special local test functions other than polynomials in Discontinuous Galerkin (DG) approaches has attracted a lot of attention and became known as DG-Trefftz methods. In particular, for the 2D Helmholtz equation plane waves have been used in \[11] to derive an Interior Penalty (IP) type Plane Wave DG (PWDG) method and to provide an a priori error analysis of its $p$-version with respect to equidistributed plane wave directions. The dependence on the distribution of the plane wave directions has been studied in \[1] based on a least squares method. However, the emphasis in \[1] has been on the $h$-version of the PWDG method, i.e., decreasing the mesh width $h$ for a fixed number $p$ of plane wave directions. In this contribution, we are interested in the $p$-version, i.e., increasing $p$ for a fixed mesh-width $h$. We formulate the choice of the plane wave directions as a control constrained optimal control problem with a continuously differentiable objective functional and the variational formulation of the PWDG method as a further constraint. The necessary optimality conditions are derived and numerically solved by a projected gradient method. Numerical results are given which illustrate the benefits of the approach.

Application of the inverse fast multipole method as a preconditioner in a 3D Helmholtz boundary element method

We investigate an efficient preconditioning of iterative methods (such as GMRES) for solving dense linear systems Ax=b that follow from a boundary element method (BEM) for the 3D Helmholtz equation, focusing on the low-frequency regime. While matrix–vector products in GMRES can be accelerated through the low-frequency fast multipole method (LFFMM), the BEM often remains computationally expensive due to the large number of GMRES iterations. We propose the application of the inverse fast multipole method (IFMM) as a preconditioner to accelerate the convergence of GMRES. The IFMM is in essence an approximate direct solver that uses a multilevel hierarchical decomposition and low-rank approximations. The proposed IFMM-based preconditioning has a tunable parameter ε that balances the cost to construct a preconditioner M, which is an approximation of A−1, and the cost to perform the iterative process by means of M. Namely, using a small (respectively, large) value of ε takes a long (respectively, short) time to construct M, while the number of iterations can be small (respectively, large). A comprehensive set of numerical examples involving various boundary value problems with complicated geometries and mixed boundary conditions is presented to validate the efficiency of the proposed method. We show that the IFMM preconditioner (with a nearly optimal ε of 10−2) clearly outperforms some common preconditioners for the BEM, achieving 1.2–10.8 times speed-up of the computations, in particular when the scale of the underlying scatterer is about five wavelengths or more. In addition, the IFMM preconditioner is capable of solving complicated problems (in a reasonable amount of time) that BD preconditioner can not.

Is Pollution Effect of Finite Difference Schemes Avoidable for Multi-Dimensional Helmholtz Equations with High Wave Numbers?

AbstractThis paper presents an approach using the method of separation of variables applied to 2D Helmholtz equations in the Cartesian coordinate. The solution is then computed by a series solutions resulted from solving a sequence of 1D problems, in which the 1D solutions are computed using pollution free difference schemes. Moreover, non-polluted numerical integration formulae are constructed to handle the integration due to the forcing term in the inhomogeneous 1D problems. Consequently, the computed solution does not suffer the pollution effect. Another attractive feature of this approach is that a direct method can be effectively applied to solve the tridiagonal matrix resulted from numerical discretization of the 1D Helmholtz equation. The method has been tested to compute 2D Helmholtz solutions simulating electromagnetic scattering from an open large cavity and rectangular waveguide.