Fast Direct Solver Research Articles

Analysis of Electromagnetic Scattering From Homogeneous Penetrable Objects by a Strong Skeletonization-Based Fast Direct Solver

This article presents a fast direct solver based on surface integral equation (SIE) to solve electromagnetic (EM) scattering from homogeneous penetrable objects. The proposed method relies on a strong admissibility skeletonization factorization (SASF) algorithm and Poggio–Miller–Chang–Harrington–Wu–Tsai (PMCHWT) formulation. In the SASF scheme, only well-separated groups that satisfy the strong admissibility condition are considered to be compressed and thus relatively fewer skeleton basis functions are selected. It is an effective way to compress the matrix with small ranks. An independent compression technique is developed for dielectric problems involving electric and magnetic currents, in which skeleton basis functions representing electric and magnetic currents are obtained separately. Moreover, a novel strategy based on matrix normalization is proposed to treat the arising “fill-in” blocks when far-field interactions are compressed. Ultimately, the impedance matrix can be cast into products of a series of block unit triangular matrices and a block diagonal matrix. Several numerical results show that the proposed approach is effective and stable as well as provides an accurate solution for scattering problems of homogeneous penetrable objects.

Fast Analysis and Optimization of Sparsely Distributed Partial Modification Problems

This article addresses the efficient numerical analysis of sparsely distributed small modifications in a large structure, that is, the sparsely distributed partial modification problem (SDPMP). A hierarchical matrix (H-matrix)-based fast direct solver is adopted to solve the SDPMP in a very short time for each new set of modifications. The main idea of this method is to reuse the parts of the compressed inverse matrix operator that do not change and only recompute the modified parts. In contrast with already existing algorithms designed to tackle a single localized modification, this new approach is still very fast for sparsely distributed small modifications, so it becomes a very attractive option to solve optimization problems efficiently, due to the very cheap analysis of the modified structure at each iteration of the optimization procedure. The efficiency of the solution of SDPMP is demonstrated for several cases involving passive and active metasurfaces and a reconfigurable antenna. As a practical example, a programmable metasurface reflector is optimized with particle swarm optimization to obtain a prescribed reflection pattern.

Preconditioned Legendre Spectral Galerkin Methods for the Non-separable Elliptic Equation

The Legendre spectral Galerkin method of self-adjoint second order elliptic equations usually results in a linear system with a dense and ill-conditioned coefficient matrix. In this paper, the linear system is solved by a preconditioned conjugate gradient (PCG) method where the preconditioner M is constructed by approximating the variable coefficients with a (T+1)-term Legendre series in each direction to a desired accuracy. A feature of the proposed PCG method is that the iteration step increases slightly with the size of the resulting matrix when reaching a certain approximation accuracy. The efficiency of the method lies in that the system with the preconditioner M is approximately solved by a one-step method based on the incomplete LU factorization technique with no fill-in, denoted by ILU(0). The ILU(0) factorization of \(M\in {\mathbb {R}}^{(N-1)^d\times (N-1)^d}\) can be computed using \({\mathcal {O}}(T^{2d} N^d)\) operations, and the number of nonzeros in the factorization factors is of \({\mathcal {O}}(T^{d} N^d)\), \(d=1,2,3\). A conclusion of the algorithm is to fast solve the resulting system from the Legendre Galerkin spectral method for Poisson equations with Dirichlet boundary conditions, which has a complexity of \({\mathcal {O}}(N^d)\). To further speed up the PCG method, an algorithm is developed for fast matrix-vector multiplications by the resulting matrix of Legendre-Galerkin spectral discretization, without the need to explicitly form it. The complexity of the fast matrix-vector multiplications is of \({\mathcal {O}}(N^d (\log _2 N)^{2})\). In view that T is independent of N in one dimension and is set to be of order \({\mathcal {O}}(\log _2 N)\) in two and three dimensions, the PCG method has a \({\mathcal {O}}(N^d (\log _2 N)^{2d})\) quasi-optimal complexity for a d dimensional domain with \((N-1)^d\) unknows, \(d=1,2,3\). In addition, a fast direct solver for the three-dimensional Poisson equation is developed, which is of \({\mathcal {O}}(N^{3} (\log _2 N)^{2})\) and improves the existing results on the computational complexity. Numerical examples are given to demonstrate the efficiency of proposed preconditioners and the algorithm for fast matrix-vector multiplications.

A New Directional Algebraic Fast Multipole Method Based Iterative Solver for the Lippmann-Schwinger Equation Accelerated with HODLR Preconditioner

We present a fast iterative solver for scattering problems in 2D, where a penetrable object with compact support is considered. By representing the scattered field as a volume potential in terms of the Green's function, we arrive at the Lippmann-Schwinger equation in integral form, which is then discretized using an appropriate quadrature technique. The discretized linear system is then solved using an iterative solver accelerated by Directional Algebraic Fast Multipole Method (DAFMM). The DAFMM presented here relies on the directional admissibility condition of the 2D Helmholtz kernel. And the construction of low-rank factorizations of the appropriate low-rank matrix sub-blocks is based on our new Nested Cross Approximation (NCA)~\cite{ arXiv:2203.14832 [math.NA]}. The advantage of our new NCA is that the search space of so-called far-field pivots is smaller than that of the existing NCAs. Another significant contribution of this work is the use of HODLR based direct solver as a preconditioner to further accelerate the iterative solver. In one of our numerical experiments, the iterative solver does not converge without a preconditioner. We show that the HODLR preconditioner is capable of solving problems that the iterative solver can not. Another noteworthy contribution of this article is that we perform a comparative study of the HODLR based fast direct solver, DAFMM based fast iterative solver, and HODLR preconditioned DAFMM based fast iterative solver for the discretized Lippmann-Schwinger problem. To the best of our knowledge, this work is one of the first to provide a systematic study and comparison of these different solvers for various problem sizes and contrast functions. In the spirit of reproducible computational science, the implementation of the algorithms developed in this article is made available at \url{https://github.com/vaishna77/Lippmann_Schwinger_Solver}.

Direct solution of the volume integral equation for electromagnetic scattering from dielectric objects with negative permittivity.

In this paper, electromagnetic scattering from dielectric objects with negative permittivity is solved by the volume integral equation (VIE)-hierarchical matrix (H-matrix)-based fast direct solver. To improve solution efficiency, a coarsening algorithm for the H-matrix is used, and the key parameter for the admissible condition is optimized through numerical experiments. Finally, the discontinuous Galerkin (DG) VIE (DGVIE) method is developed for the fast direct solver. It is shown that scattering from inhomogeneous dielectric objects with negative permittivity and with multiscale structure can be analyzed efficiently by the DGVIE-H-matrix-based fast direct solver.

Zeta correction: a new approach to constructing corrected trapezoidal quadrature rules for singular integral operators

A high-order accurate quadrature rule for the discretization of boundary integral equations (BIEs) on closed smooth contours in the plane is introduced. This quadrature can be viewed as a hybrid of the spectral quadrature of Kress (Math. Comput. Model. 15(3-5), 229–243 1991) and the locally corrected trapezoidal quadrature of Kapur and Rokhlin (SIAM J. Numer. Anal. 34(4), 1331–1356, 1997). The new technique combines the strengths of both methods, and attains high-order convergence, numerical stability, ease of implementation, and compatibility with the “fast” algorithms (such as the Fast Multipole Method or Fast Direct Solvers). Important connections between the punctured trapezoidal rule and the Riemann zeta function are introduced, which enable a complete convergence analysis and lead to remarkably simple procedures for constructing the quadrature corrections. The paper reports a detailed comparison between the new method and the methods of Kress, of Kapur and Rokhlin, and of Alpert (SIAM J. Sci. Comput. 20(5), 1551–1584, 1999).

Strong Admissibility Skeletonization Factorization for Fast Direct Solution of Electromagnetic Scattering From Conducting Objects

In this article, a fast direct solver based on strong admissibility skeletonization factorization (SASF) is proposed for electromagnetic scattering from conducting objects. Different from the conventional skeletonization scheme, the proposed skeletonization constructs the hierarchical matrix representation, in which only far-field interactions are compressed. As a result, the approximation rank is relatively small, and the computational efficiency will have significant improvement. Subsequently, the strong skeletonization factorization is applied to the compressed system matrix. The system matrix can be factorized into products of a series of block unit triangular matrices and a block diagonal matrix. The arising fill-in blocks corresponding to far-field interactions are compressed and eliminated by a novel and efficient method to maintain the high efficiency and accuracy of the factorization procedure. The computational complexity and storage requirement of the proposed factorization scale as <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$O(N^{1.5})$ </tex-math></inline-formula> and <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$O(N \log N)$ </tex-math></inline-formula> , respectively. Several numerical results are presented to demonstrate the accuracy and effectiveness of the proposed method.

Fast direct solver for CN-ADI-FV scheme to two-dimensional Riesz space-fractional diffusion equations

Based on the Crank-Nicolson alternating direction implicit finite volume (CN-ADI-FV) discretization of two dimensional Riesz space-fractional diffusion equations (2 D RSFDEs) proposed in [1], a fast direct solver is proposed to obtain the numerical solutions. Linear systems result from the CN-ADI-FV discretization can be treated as two matrix equations with symmetric positive definite (SPD) Toeplitz coefficient matrices. Then numerical solutions, expressed as multiple matrix-vector multiplications involving the inverse of SPD Toeplitz matrices, can be effectively obtained by using Gohberg-Semencul formula (GSF), where the preconditioned conjugate gradient (PCG) method is applied as an acceleration technique. Moreover, under certain assumptions, the spectral of the preconditioned matrices of the linear system associated with the GSF are theoretically proved to be clustered around 1. Then the PCG method is verified to converge superlinearly and can be implemented with O(NlogN) operations, where N is the number of spatial unknowns in both x- and y- directions. Besides, the PCG method is used only twice in the whole solution process, and the complete solution process can be finished with O(MN2logN) operations and O(N) memory, where M is the number of temporal unknowns. Numerical experiments are given to demonstrate the efficiency of our approach.

An alternative extended linear system for boundary value problems on locally perturbed geometries

This manuscript presents a new extended linear system for integral equation based techniques for solving boundary value problems on locally perturbed geometries. The new extended linear system is similar to a previously presented technique for which the authors have constructed a fast direct solver. The key features of the work presented in this paper are that the fast direct solver is more efficient for the new extended linear system and that problems involving specialized quadrature for weakly singular kernels can be easily handled. Numerical results illustrate the improved performance of the fast direct solver for the new extended system when compared to the fast direct solver for the original extended system.

A Low-Rank Schwarz Method for Radiative Transfer Equation With Heterogeneous Scattering Coefficient

Random sampling has been used to find low-rank structure and to build fast direct solvers for multiscale partial differential equations of various types. In this work, we design an accelerated Schw...

Interior eigensolver for sparse Hermitian definite matrices based on Zolotarev’s functions

This paper proposes an efficient method for computing selected generalized eigenpairs of a sparse Hermitian definite matrix pencil $(A,B)$. Based on Zolotarev's best rational function approximations of the signum function and conformal mapping techniques, we construct the best rational function approximation of a rectangular function supported on an arbitrary interval via function compositions with partial fraction representations. This new best rational function approximation can be applied to construct spectrum filters of $(A,B)$ with a smaller number of poles than a direct construction without function compositions. Combining fast direct solvers and the shift-invariant generalized minimal residual method, a hybrid fast algorithm is proposed to apply spectral filters efficiently. Compared to the state-of-the-art algorithm FEAST, the proposed rational function approximation is more efficient when sparse matrix factorizations are required to solve multi-shift linear systems in the eigensolver, since the smaller number of matrix factorizations is needed in our method. The efficiency and stability of the proposed method are demonstrated by numerical examples from computational chemistry.

A fast direct solver for nonlocal operators in wavelet coordinates

In this article, we consider fast direct solvers for nonlocal operators. The pivotal idea is to combine a wavelet representation of the system matrix, yielding a quasi-sparse matrix, with the nested dissection ordering scheme. The latter drastically reduces the fill-in during the factorization of the system matrix by means of a Cholesky decomposition or an LU decomposition, respectively. This way, we end up with the exact inverse of the compressed system matrix with only a moderate increase of the number of nonzero entries in the matrix.To illustrate the efficacy of the approach, we conduct numerical experiments for different highly relevant applications of nonlocal operators: We consider (i) the direct solution of boundary integral equations in three spatial dimensions, issuing from the polarizable continuum model, (ii) a parabolic problem for the fractional Laplacian in integral form and (iii) the fast simulation of Gaussian random fields.

Accelerated Direct Solver With Multiscale Compressed Block Decomposition and Multilevel Characteristic Basis Function Method

In this letter, an accelerated fast direct solver based on the multiscale compressed block decomposition (MSCBD) and the multilevel characteristic basis function method (MLCBFM) is presented to analyze electrically large problems. The MSCBD is used to enhance the efficiency of the characteristic basis functions generation to accelerate the process of linear system solving in the MLCBFM. In addition, the MLCBFM can reduce the number of unknowns and improve the efficiency of MSCBD process. That means the MSCBD and MLCBFM can accelerate each other. The correctness and efficiency of the method are demonstrated with numerical results for several perfect electric conducting objects, by comparing with other methods.

A fast direct solver for two dimensional quasi-periodic multilayered media scattering problems

This manuscript presents a fast direct solution technique for solving two dimensional wave scattering problems from quasi-periodic multilayered structures. When the interface geometries are complex, the dominant term in the computational cost of creating the direct solver scales $O(NI)$ where $N$ is the number of discretization points on each interface and $I$ is the number of interfaces. The bulk of the precomputation can be re-used for any choice of incident wave. As a result, the direct solver can solve over 200 scattering problems involving an eleven layer geometry with complex interfaces 100 times faster than building a new fast direct solver from scratch for each new set of boundary data. An added benefit of the presented solver is that building an updated solver for a new geometry involving a replaced interface or a change in material property in one layer is inexpensive compared to building a new fast direct solver from scratch.

Hybrid FEM–SBM solver for structural vibration induced underwater acoustic radiation in shallow marine environment

This study focuses on the efficient computational model for predicting the structural vibration induced underwater acoustic radiation in shallow marine environment. In the proposed model, the finite element method (FEM) is used to calculate the vibration response of the shell structures and the singular boundary method (SBM) with near-field and far-field Green’s functions is used to simulate the underwater acoustic radiation excited by shell structural vibration in shallow water marine environment, namely, the hybrid FEM–SBM solver. For axisymmetric shell structures, the axisymmetric numerical formulation of the FEM–SBM solver is derived to save the computational cost. For general shell structures, the fast direct solver based on the recursive skeletonization factorization is introduced into the FEM–SBM solver to accelerate the computation of the underwater acoustic field. Numerical results demonstrate that the proposed FEM–SBM solver provides more efficient and reliable solutions in comparison with the coupled FEM/FEM solver (COMSOL software). The effect of the sound velocity profiles on the underwater acoustic radiation field in shallow water marine environment is numerically investigated.

Direct Solution of Scattering Problems Using Generalized Source Integral Equations

A class of inherently compressible integral equation formulations for problems of scattering by impenetrable objects, which makes use of generalized directional sources, is presented. The new formulation effectively reduces the problem's dimensionality and, thus, allows for efficient low-rank compression of moment matrices' off-diagonal blocks. When the formulation is used with a hierarchical matrix compression and factorization algorithm, a fast direct solver is obtained. The computational bottlenecks introduced by the proposed generalized formulation, in both the matrix-fill and matrix compression stages, are alleviated by using nonuniform sampling-based techniques. These techniques are described in detail for one choice of generalized sources, which use absorbing equivalent source shields, and can be extended to other shield types. The formulation's properties and limitations are studied and its enhanced compressibility is used for the development of a fast direct solver.

A novel fast direct solver for 3D elastic inclusion problems with the isogeometric boundary element method

We present a novel fast direct solver to simulate 3D large scale elastic inclusion problems. The method combines the isogeometric analysis boundary element method (IGABEM) and the hierarchical off-diagonal low-rank (HODLR) matrix based on non-uniform rational B-splines (NURBS). Hence the 3D geometric surface can be accurately described by the bivariate NURBS basis functions. In order to solve the large scale problems, a stable accelerated algorithm is used to approximate the off-diagonal submatrices by low-rank matrices. Based on the accelerated algorithm, a hybrid approximation algorithm consisting of singular value decomposition (SVD) and adaptive cross approximation (ACA) is proposed to solve the 3D elastic inclusion problems. The validity and accuracy of the method are verified by testing the four methods. Among the numerical results obtained from the four methods, the method proposed in this paper uses less CPU time and storage space to obtain accurate results.

A fast Fourier transform based direct solver for the Helmholtz problem

SummaryThis article is devoted to the efficient numerical solution of the Helmholtz equation in a two‐ or three‐dimensional (2D or 3D) rectangular domain with an absorbing boundary condition (ABC). The Helmholtz problem is discretized by standard bilinear and trilinear finite elements on an orthogonal mesh yielding a separable system of linear equations. The main key to high performance is to employ the fast Fourier transform (FFT) within a fast direct solver to solve the large separable systems. The computational complexity of the proposed FFT‐based direct solver is operations. Numerical results for both 2D and 3D problems are presented confirming the efficiency of the method discussed.

GPU acceleration of CaNS for massively-parallel direct numerical simulations of canonical fluid flows

This work presents the GPU acceleration of the open-source code CaNS for very fast massively-parallel simulations of canonical fluid flows. The distinct feature of the many-CPU Navier–Stokes solver in CaNS is its fast direct solver for the second-order finite-difference Poisson equation, based on the method of eigenfunction expansions. The solver implements all the boundary conditions valid for this type of problems in a unified framework. Here, we extend the solver for GPU-accelerated clusters using CUDA Fortran. The porting makes extensive use of CUF kernels and has been greatly simplified by the unified memory feature of CUDA Fortran, which handles the data migration between host (CPU) and device (GPU) without defining new arrays in the source code. The overall implementation has been validated against benchmark data for turbulent channel flow and its performance assessed on a NVIDIA DGX-2 system (16 T V100 32Gb, connected with NVLink via NVSwitch). The wall-clock time per time step of the GPU-accelerated implementation is impressively small when compared to its CPU implementation on state-of-the-art many-CPU clusters, as long as the domain partitioning is sufficiently small that the data resides mostly on the GPUs. The implementation has been made freely available and open source under the terms of an MIT license.

An efficient Navier-Stokes based numerical wave tank using fast Poisson solvers and the immersed boundary method

An efficient numerical model for wave induced flows and their interaction with impermeable coastal structures, is presented. Using recent developments in the field of two-fluid (water/air) Navier-Stokes (NS) simulations, an efficient numerical wave tank (NWT) was developed which combines the Immersed Boundary (IB) method and Fast Direct Solvers (FDS). Conventional NS-based NWTs require significant computational resources which prevents their use in discretisation demanding cases or extensive parametric studies. The computational efficiency of the new NWT was achieved by replacing the conventional variable-coefficient Poisson equation with a constant-coefficient one which is solved with a FDS. Results of the accuracy, limitations, and computational performance of the NWT are presented and discussed based on several validation tests: (1) propagation of a 2nd-order Stokes wave, (2) formation of a standing wave, (3) spilling breaking over a mild slope bed, (4) three-dimensional (3D) interaction of a solitary wave with a vertical abutment, and (5) diffraction of a solitary wave by a single-row pile breakwater. The computational performance of the new NWT was demonstrated for problems with 1–280 million grid nodes. The overall speed-up of the new NWT increases with the problem size, and it is estimated to be about 30 for large 3D problems.