A Linear Complexity Direct Volume Integral Equation Solver for Full-Wave 3-D Circuit Extraction in Inhomogeneous Materials

We develop a direct finite-element solver of linear (optimal) complexity to extract broadband circuit parameters such as S-parameters of arbitrarily shaped 3-D interconnects in inhomogeneous dielectrics. Numerical experiments demonstrate a clear advantage of the proposed solver as compared with existing finite-element solvers that employ state-of-the-art direct sparse matrix solutions. A linear complexity in both CPU time and memory consumption is achieved with prescribed accuracy satisfied. A finite-element matrix from the analysis of a large-scale 3-D circuit in multiple dielectrics having 5.643 million unknowns is directly factorized in less than 2 hours on a single core running at 2.8 GHz.

In this article, the 3-D electromagnetic (EM) scattering by 2-D dielectric arbitrary anisotropic scatterers embedded in a layered uniaxial anisotropic medium is studied. This 2.5-D EM scattering problem is mathematically formulated by the volume integral equation (VIE) whose discretized weak forms are based on the 2.5-D roof-top basis functions and solved by the stabilized biconjugate gradient fast Fourier transform (BCGS-FFT). Meanwhile, the 2.5-D dyadic Green’s functions (DGFs) for the layered uniaxial media are given in detail and their evaluation is also discussed. In particular, a tricky variable replacement strategy is proposed to obtain analytical expressions for partial 2.5-D DGF components. Besides the validation of the 2.5-D DGFs by comparing them with the corresponding 3-D values, several numerical experiments are also carried out to validate the accuracy and efficiency of the BCGS-FFT solver for the 2.5-D EM scattering in the layered anisotropic circumstance by comparing the results with those obtained by a 3-D VIE solver. The major new contribution of this work is to extend the 2.5-D EM scattering computation to accommodate the uniaxial anisotropy of the layered background medium and the arbitrary anisotropic scatterers.

Summary form only given. Fast direct integral equations methods for analyzing time-harmonic scattering phenomena constitute an active area of research in the CEM community. These solvers present an interesting alternative to FMM-based iterative schemes for problems that are inherently ill-conditioned and/or involve many excitations. Present fast direct methods leverage the low-rank (LR) nature of blocks of the LU factors or inverse of the method of moment (MoM) matrix (if and when this property manifests itself). In essence, they leverage the limited number of degrees of freedom of the fields radiated by sources in the presence of the scatterer (or portions thereof), observed across a finite region of the scatterer (whenever the concept makes sense). LR compression strategies do not lead to low-complexity solvers for electrically large concave structures. This fact was demonstrated in [H. Guo et al., Antennas and Wireless Propagation Letters, IEEE, vol. 11, 2012], where blocks of LU factors that proved incompressible by LR methods were recompressed by MLMDA/butterfly methods. Somewhat unexpectedly, this study showed that MLMDA/butterfly scheme produced compression above and beyond that provided by LR schemes, requiring only O(N log2 N) final storage irrespective of the nature or electrical size of the scatterer. The practical applicability of the scheme however was severely hindered by the fact that the MLMDA/butterfly scheme acted as an add-on to a LR solver. Here, we report on a new direct integral equation solver by that utilizes MLMDA/butterfly schemes for compressing the LU factors of a MoM matrix. When compared to its predecessor, the new solver has three important features of note. (i) The solver entirely bypasses the LR compression step of its predecessor and operates directly on butterfly-compressed blocks. The latter is achieved using new schemes for rapidly adding and multiplying butterfly-compressed operators. (ii) The solver executes in parallel, using a hybrid OpenMP-MPI strategy to accelerate the hierarchical inverse/decomposition of the MoM matrix. Despite the inherently sequential nature of this operation, the solver exhibits excellent scaling properties up to several hundred processors. (iii) The solver applies to large-scale 3D (as opposed to 2D) analysis and was implemented to invert a combined field (as opposed to an electric field) integral operator. The solver is capable of inverting MoM matrices that discretize combined field integral equations modeling scattering from electrically large structures involving millions of unknowns on a small cluster using just a few hours of CPU time.

Summary form only given. Solutions of large-scale electromagnetic scattering problems are typically obtained with Rao-Wilton-Glisson (RWG) Method of Moments (MoM) accelerated with Multi-Level-Fast-Multipole-Algorithm (MLFMA). Due to the low-order RWG MoM discretization and iterative nature of MLFMA such solutions present various challenges when realistic large-scale targets are considered. The most salient of these challenges are: 1) inefficient error control of the solution beyond 2-3 digits, 2) poor conditioning of the matrix resulting from multiscale discretization of the model and/or the presence of highly resonant phenomena, and 3) the necessity to repeat the solution process for each excitation. In this work, we present a novel computational framework, which efficiently overcomes these three challenges. Our method is based on Non-Uniform-Rational-B-Splines (NURBS) representation of the scatterer geometries and Higher-Order (HO) Locally Corrected Nystrom (LCN) discretization of the pertinent Combined Field Integral Equation (CFIE). Due to the HO nature of the LCN scheme, scattering solutions of desired precision are obtained with exponentially higher efficiency as compared to the RWG MoM (Jeffrey, et. al., IEEE AP Mag. pp. 294-308, vol. 55, no. 3, June 2013). This allows us to accurately reconstruct the fields in both lit and deeply shaded areas, eliminate lack of causality in time domain responses synthesized from broadband frequency domain solutions, and overcome various other problems stemming from insufficient solution accuracy. The second and third challenges are addressed through fast direct solution of the matrix equation resultant from LCN discretization of the CFIE. Our fast direct solution of the matrix equation is based on the block-LU decomposition aided with Adaptive-Cross-Approximation (ACA) compression of the blocks similar to Shaeffer's algorithm (Shaeffer, IEEE TAP, no. 8, pp. 2306-2313, 2008). To overcome the computational complexity of the fast direct solution, the block LU decomposition is performed on multiple Graphics Processing Units (GPUs) with fast Out-of-Core memory augmentation. A commodity workstation featuring two Xeon hexacore processors, 2 K10 GPUs, 256Gb of RAM, and SSD RAID array, affords us solutions for a broad range of complex scattering problems, previously demonstrated with over 1 million RWG degrees of freedom, which cannot be obtained using iterative algorithms. In the case of problems requiring solutions with multiple excitations such as antenna array modeling and scattering problems with many incidence angles, despite the large computational complexity of the direct block LU factorization, the overall solution time was seen to be substantially lower than in the case of a MLFMA accelerated iterative approach as the latter must be repeated for each excitation.

The fast direct solver based on skeletonization is extended to the method of moments (MoM) solution of volume integral equations (VIE). The impedance matrix is approximated by the hierarchical data sparse representation through skeletonization. The inverse matrix is then implicitly computed. Numerical experiments are conducted to study the performance of the direct solver.

A novel double‐higher‐order entire‐domain volume integral equation (VIE) technique for efficient analysis of electromagnetic structures with continuously inhomogeneous dielectric materials is presented. The technique takes advantage of large curved hexahedral discretization elements—enabled by double‐higher‐order modeling (higher‐order modeling of both the geometry and the current)—in applications involving highly inhomogeneous dielectric bodies. Lagrange‐type modeling of an arbitrary continuous variation of the equivalent complex permittivity of the dielectric throughout each VIE geometrical element is implemented, in place of piecewise homogeneous approximate models of the inhomogeneous structures. The technique combines the features of the previous double‐higher‐order piecewise homogeneous VIE method and continuously inhomogeneous finite element method (FEM). This appears to be the first implementation and demonstration of a VIE method with double‐higher‐order discretization elements and conformal modeling of inhomogeneous dielectric materials embedded within elements that are also higher (arbitrary) order (with arbitrary material‐representation orders within each curved and large VIE element). The new technique is validated and evaluated by comparisons with a continuously inhomogeneous double‐higher‐order FEM technique, a piecewise homogeneous version of the double‐higher‐order VIE technique, and a commercial piecewise homogeneous FEM code. The examples include two real‐world applications involving continuously inhomogeneous permittivity profiles: scattering from an egg‐shaped melting hailstone and near‐field analysis of a Luneburg lens, illuminated by a corrugated horn antenna. The results show that the new technique is more efficient and ensures considerable reductions in the number of unknowns and computational time when compared to the three alternative approaches.

Parallel direct method for linear equations

Time-domain layered finite element reduction recovery (LAFE-RR) method was recently developed for large-scale electromagnetic analysis of high-speed integrated circuits (ICs). This method is capable of analytically and rigorously reducing the system matrix of a 3-D multilayer circuit to that of a single-layer one regardless of the original problem size. In addition, the reduced system matrix preserves the sparsity of the original system matrix. In this paper, an efficient algorithm is proposed to recover the volume unknowns in the time-domain LAFE-RR method. This algorithm constitutes a direct solution of the matrix formed by volume unknowns in each layer. This direct solution possesses a linear complexity in both central processing unit (CPU) time and memory consumption. The cost of matrix inversion is negligible. The cost of matrix solution scales linearly with the matrix size. Numerical and experimental results have demonstrated the accuracy and efficiency of the proposed algorithm.

A butterfly-based direct combined-field integral equation (CFIE) solver for analyzing scattering from electrically large, perfect electrically conducting objects is presented. The proposed solver leverages the butterfly scheme to compress blocks of the hierarchical LU-factorized discretized CFIE operator and uses randomized butterfly reconstruction schemes to expedite the factorization. The memory requirements and computational cost of the direct butterfly-CFIE solver scale as $O(N\mathrm{log}^2N)$ and $O(N^{1.5}\mathrm{log}N)$, respectively. These scaling estimates permit significant memory and CPU savings when compared to those realized by low-rank (LR) decomposition-based solvers. The efficacy and accuracy of the proposed solver are demonstrated through its application to the analysis of scattering from canonical and realistic objects involving up to 14 million unknowns.

In this work, we prove that the sparse matrix resulting from a finite-element-based analysis of electrodynamic problems can be represented by an <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">H</i> matrix without any approximation, and the inverse of this sparse matrix has a data-sparse <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">H</i> -matrix approximation with error well controlled. Two proofs are developed. One is based on the general eigenvalue-based solution to the ordinary differential equations, and the other is based on the relationship between a partial differential operator and an integral operator. Both proofs have reached the same conclusion. Based on the proof, we develop an <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">H</i> -matrix-based direct finite-element solver of <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">O</i> ( <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">kN</i> log <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">N</i> ) memory complexity and <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">O</i> ( <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">k</i> <sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">2</sup> <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">N</i> log <sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">2</sup> <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">N</i> ) time complexity for solving electromagnetic problems, where <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">k</i> is a small parameter that is adaptively determined based on accuracy requirements, and <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">N</i> is the number of unknowns. Both inverse-based and LU-based direct solutions are developed. The LU-based solution is further accelerated by nested dissection. A comparison with the state-of-the-art direct finite element solver that employs the most advanced sparse matrix solution has shown clear advantages of the proposed direct solver. In addition, the proposed solver is applicable to arbitrarily-shaped three-dimensional structures and arbitrary inhomogeneity.

A direct domain-decomposition-based time-domain finite-element method of linear complexity is developed to overcome the challenge of simulating a wide range of geometrical scales present in an integrated circuit (IC) system. Via a set of orthogonal prism vector basis functions, we decompose the entire 3-D system of an IC problem into 2-D subsystems and then into 1-D subsystems with negligible computational cost. We then further decompose each 1-D subsystem into two domains, one with conductors and the other one with dielectrics. A direct solution without any iteration is proposed to solve the resulting system in linear complexity. Furthermore, the method allows for the use of different meshes in different domains. Numerical simulations of ICs and package problems have demonstrated the accuracy, linear complexity, and meshing flexibility of the proposed method.

Multi-level numerical methods that obtain the exact solution of a linear system are presented. The methods are devised by combining ideas from the full multi-grid algorithm and perfect reconstruction filters. The problem is stated as whether a direct solver is possible in a full multi-grid scheme by avoiding smoothing iterations and using different coarse grids at each step. The coarse grids must form a partition of the fine grid and thus establishes a strong connection with domain decomposition methods. An important analogy is established between the conditions for direct solution in multi-grid solvers and perfect reconstruction in filter banks. Furthermore, simple solutions of these conditions for direct multi-grid solvers are found by using mirror filters. As a result, different configurations of direct multi-grid solvers are obtained and studied.

A linear-complexity direct matrix solution is developed for the surface-integral based impedance extraction of arbitrarily-shaped 3-D non-ideal conductors embedded in dielectric materials. It outperforms state-of-the-art impedance solvers with fast CPU-time, modest memory-consumption, and without sacrificing accuracy. The inverse of a 2.6-million-unknown matrix arising from the extraction of large-scale 3-D interconnects was obtained in 1.5 GB memory and 1.3 hours on a 3 GHz CPU.

This paper discusses compressive sensing with sparse measurement matrices. Sparse matrices have several attractive properties, like low computational complexity in both encoding and recovery, easy incremental updates to signals, and low storage requirement, etc. Typical examples of existing algorithms for sparse signal recovery with sparse measurement matrices include convex relaxation, matching pursuit, and Bayesian framework based approaches. In this paper, we propose an alternative technique to this problem. The proposed technique has a linear recovery complexity and relatively good empirical behavior. In this technique, we employ a permutation-based multi-dimensional measurement matrix, which is composed of several sub-matrices, each consisting of a block-diagonal matrix and a random permutation matrix. The measurement symbols generated from the same permutation matrix are referred to as a dimension. Such a measurement matrix brings some useful features to the measurement symbols. Fully exploiting these features enables us to design a simple and effective recovery algorithm. The proposed recovery algorithm employs an iterative process. In each iteration, the algorithm looks for measurement symbols with certain features, and uses such features to recover the source symbols related to these measurements. An iterative process is then applied, along with an interference cancellation operation, to reconstruct all source symbols gradually. The complexity of the proposed algorithm is relatively low, which grows linearly with the source signal length. Numerical results show that the proposed technique empirically offers a much lower sketch length than ℓ <sub xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">1</sub> -minimization-based convex relaxation and Bayesian framework based algorithms. It achieves the empirical lower bound of sketch length and linear recovery complexity at the same time.

${\mathcal{ H}}^{2}$ -matrix constitutes a general mathematical framework for efficient computation of both partial-differential-equation (PDE) and integral-equation (IE)-based operators. Existing linear-complexity ${\mathcal{ H}}^{2}$ matrix-matrix product (MMP) algorithm lacks explicit accuracy control, while controlling accuracy without compromising linear complexity is challenging. In this article, we develop an accuracy controlled ${\mathcal{ H}}^{2}$ MMP algorithm by instantaneously changing the cluster bases during the matrix product computation based on prescribed accuracy. Meanwhile, we retain the computational complexity of the overall algorithm to be linear. Different from the existing ${\mathcal{ H}}^{2}$ MMP algorithm where formatted multiplications are performed using the original cluster bases, in the proposed algorithm, all additions and multiplications are either exact or computed based on prescribed accuracy. Furthermore, the original ${\mathcal{ H}}^{2}$ -matrix structure is preserved in the matrix product. While achieving optimal complexity for constant-rank matrices, the computational complexity of the proposed algorithm is also minimized for variable-rank ${\mathcal{ H}}^{2}$ -matrices. For example, it has a complexity of $O(NlogN)$ for computing electrically large volume IEs, where $N$ is matrix size. The proposed work serves as a fundamental arithmetic in the development of fast solvers for large-scale electromagnetic analysis. Applications to both large-scale capacitance extraction and electromagnetic scattering problems involving millions of unknowns on a single core have demonstrated the accuracy and efficiency of the proposed algorithm.