1 Introduction↩︎

The numerical solution of large-scale electromagnetic scattering problems is important in both physics and engineering applications. For problems defined over bounded domains, numerical techniques such as the finite difference method and the finite element method are effective. However, to solve problems in unbounded domains, these methods become challenging due to the necessity of handling radiation conditions. Consequently, numerical solvers based on boundary integral equations are a promising alternative.

The fast multipole method (FMM) is well known as an efficient acceleration technique of boundary integral solvers for scattering problems [1]. Parallelization techniques for FMM in distributed memory systems have also been well studied, for example, [2]–[4]. The FMM is typically interpreted as a method for speeding up matrix-vector multiplication, and is therefore often coupled with Krylov subspace iterative solvers like GMRES. Consequently, when solving problems involving many right-hand sides, the required computational time, in the worst-case scenario, scales as a multiple of the number of right-hand sides. The inherent drawback of using iterative solvers with FMM can be effectively handled by employing fast direct solvers (FDS) [5].

This study focuses on the transverse electric (TE)-mode electromagnetic scattering in transmission problems. Because TE-mode scattering reduces to two-dimensional Helmholtz scattering, using solvers based on the hierarchically semiseparable (HSS) representation [6] is attractive. FDSs based on the HSS representation include the Martinsson–Rokhlin solver [7], HSS-ULV factorization [8], and recursive skeletonization [9]. These HSS solvers exhibit high parallel efficiency because all off-diagonal blocks are recursively approximated as low rank, although alternative representations, such as strongly admissible recursive skeletonization [10] and \(\mathcal{H}^2\)-ULV factorization [11], are required for three-dimensional problems. However, the HSS representation is not suitable in high-order discretization because local corrections or a support of basis functions inevitably spill over into adjacent off-diagonal blocks even when solving two-dimensional problems like that shown in Figure 1, leading to a loss of accuracy of low-rank approximation. Therefore, a major limitation for HSS solvers remains: we are forced to compromise on the discretization order to retain high parallel efficiency.

This study developed an MPI/OpenMP parallelized scalable HSS-type FDS with \(O(h^3)\) convergence for TE-mode transmission problems, whereas conventional HSS solvers typically yield only \(O(h)\) convergence (where \(h\) denotes the quadrature node spacing). This work makes two main contributions: a boundary integral formulation tailored to the HSS representation and an efficient load-balancing technique.

The developed HSS solver is based on the weakly singular Burton–Miller (BM) boundary integral equation, which is expected to avoid quasi-resonances [12]. Due to this mild singularity, the solver achieves \(O(h^3)\) convergence via the Nyström method with a one-point corrected quadrature, in which only the diagonal terms are locally corrected. Although an efficient parallelization technique for the HSS-ULV factorization has been described [13], we observed that the dominant computational bottleneck lies in calculating matrix entries rather than the factorization itself. These entries are computed during low-rank approximations, for which we employ the proxy-surface method [7]. While the proxy-surface method has been successfully extended to the Helmholtz [14], Maxwell’s [15], and Navier–Cauchy [16] equations, its parallel load balancing in an HSS solver requires further optimization. Therefore, we use an appropriate block formulation for the linear system to balance the computational load across processes and threads. Numerical results demonstrate that our HSS direct solver achieves high-order convergence and nearly ideal strong and weak scalabilities.

2 Boundary integral equations for TE-mode transmission scattering problems↩︎

Let us consider the Cartesian coordinate system \((x_1, x_2, x_3)\). There is a bounded scatterer made of a penetrable object with a smooth surface. This scatterer is shaped as a general cylinder that extends infinitely along the \(x_3\) direction. A time harmonic electromagnetic wave propagates along the \(x_3\) axis. When the electric field component is entirely perpendicular to the direction of wave propagation, this specific propagation is called the TE-mode. Under these assumptions, the scattering problem is reduced to the transmission problem of the Helmholtz equation in the \(x_1\)–\(x_2\) plane. Let \(\Omega_1 \subset \mathbb{R}^2\) be the part of the scatterer cut off at this \(x_1\)–\(x_2\) plane. We define \(\Omega_0 = \mathbb{R}^2 \setminus \Omega_1\) as the exterior connected domain. For a given function \(v\), \(S_{k}v\) and \(D_{k}v\) are the single- and double-layer potentials with wavenumber \(k\), and \(D_{k}^* v\) and \(N_{k} v\) are their normal derivatives, respectively. To solve the TE-mode scattering, the following boundary integral equation (BIE) can be used: \[\require{physics} \begin{gather} \mqty( -\qty(\varepsilon_{0} S_{k_{0}} + \varepsilon_{1} S_{k_1}) & D_{k_0} + D_{k_1} \\ -\qty(D_{ k_0}^{\ast} + D_{k_1}^{\ast}) & \frac{1}{\varepsilon_{0}} N_{k_0} + \frac{1}{\varepsilon_{1}} N_{k_1} ) \mqty( q(x) \\ u(x) ) \\ = \mqty( -u^{\mathrm{in}}(x) \\ -\frac{1}{\varepsilon_{0}} \pdv{u^{\mathrm{in}}}{n}\qty(x) ), \quad x \in \Gamma, \label{eq:pmchwt} \end{gather}\tag{1}\] where \(\Gamma = \partial \Omega_1\) is the boundary of a general cylindrical domain in the \(x_1\)–\(x_2\) plane, \(n(z)\) is the outward unit normal vector at \(z \in \Gamma\), \(u\) is the limit to \(\Gamma\) of the magnetic field in the \(x_3\) direction, \(u^{\mathrm{in}}\) is the incident wave, and \(q\) is the scaled normal derivative of \(u\). For \(j=0, 1\), the wavenumber \(k_j\) in \(\Omega_j\) is determined by \(k_j = \omega \sqrt{\varepsilon_j \mu_j}\), where \(\omega\) is the angular frequency of the incident wave and \(\varepsilon_j\) and \(\mu_j\) are the constant permittivity and permeability, respectively. The above BIE, known as the PMCHWT equation, is the standard formulation for transmission problems [1], and it can be discretized using the one-point zeta-corrected quadrature [17].

As an alternative formulation for the same transmission problem, the BM equation is given as \[\require{physics} \begin{gather} \mqty( D_{k_1} + I/2 & -\varepsilon_1 S_{k_1} \\ D_{k_0} - I/2 + \alpha N_{k_0} & -\varepsilon_0 \qty{S_{k_0} +\alpha \qty(D_{k_0}^* + I/2)} ) \\ \cdot \mqty( u(x) \\ q(x) ) = \mqty( 0 \\ -u^{\mathrm{in}}(x) - \alpha \pdv{u^{\mathrm{in}}}{n}\qty(x) ), \quad x \in \Gamma, \label{eq:bm} \end{gather}\tag{2}\] where \(I\) is the identity operator and \(\alpha = i/k_0\) is a constant of the BM method. While the PMCHWT equation 1 involves eight layer potentials, the BM equation contains only six. This reduction in the number of layer potentials leads to an efficient FDS [18]. However, a solver based on these two equations achieves only \(O(h)\) convergence in the one-point zeta-corrected quadrature, because \(N_k v\) possesses a hypersingularity, whereas the other layer potentials \(S_k v\), \(D_k v\), and \(D_k^* v\) are only weakly singular.

We therefore employ the weakly singular BM equation, defined as \[\require{physics} \begin{align} \left( \begin{array}{c:c} D_{k_1} + I/2 & -\varepsilon_1 S_{k_1} \\ \hdashlineD_{k_0} &- I/2 \\ &+ \alpha \qty(N_{k_0} - N_{k_1})& \begin{align} -\epsilon_0 S_{k_0} &- \alpha \qty(\epsilon_0 D_{k_0}^* - \epsilon_1 D_{k_1}^* ) \\ &- \alpha (\varepsilon_0 + \varepsilon_1)I/2 \end{align} \end{array} \right) \\ \cdot \mqty( u(x) \\ q(x) ) = \mqty( 0 \\ -u^{\mathrm{in}}(x) - \alpha \pdv{u^{\mathrm{in}}}{n}\qty(x) ), \quad x \in \Gamma. \label{eq:weakly95bm} \end{align}\tag{3}\] This formulation is equivalent to the one described in [19]. Because the difference \((N_{k_0} - N_{k_1})v\) has at most a weakly singular kernel, this achieves \(O(h^3)\) convergence even with the one-point corrected quadrature. We note that, once the boundary values \(u\) and \(q\) are obtained, the solution inside the domain can be calculated through integral representations.

3 Distributed parallel fast direct solver↩︎

This study developed a load-balanced, MPI/OpenMP hybrid parallel FDS. This solver is a variant of the Martinsson–Rokhlin solver [16] and features two key characteristics that improve load balancing:

• The adoption of coarse-grained MPI processes to minimize the volume of MPI communications. Unlike the FMM, which involves the communication of vectors only, this solver requires the exchange of both vectors and matrices between processes.

• The introduction of a block formulation to ensure that all OpenMP threads handle the same number of layer potentials. As shown in 3 , the number of layer potentials is non-uniform across the \(2 \times 2\) blocks.

We note that the direct solver follows the standard procedure except for its parallelization scheme.

An overview of the parallelization strategy is shown in Figure 2. First, the boundary \(\Gamma\) is recursively partitioned using a perfect binary tree structure. The entire tree is partitioned into subtrees. According to this partitioning, the entire index \(I_0\) is partitioned into \(p\) subsets \(I_1, I_2, \ldots, I_p\). Each MPI process owns a single subtree. The linear system arising from 1 , 2 , and 3 can be discretized into the following block forms via partitioning of the index: \[\mqty( A_{11} & A_{12} & \cdots & A_{1p} \\ A_{21} & A_{22} & \cdots & A_{2p} \\ \vdots & \vdots &\ddots &\vdots \\ A_{p1} & A_{p2} & \cdots & A_{pp} ) \mqty( \varphi_1 \\ \varphi_2 \\ \vdots \\ \varphi_p ) = \mqty( b_1 \\ b_2 \\ \vdots \\ b_p ), \label{eq:block95linear}\tag{4}\] where \(\varphi_j\) and \(b_j\) correspond to the discretized solution and right-hand side, respectively. Each block \(A_{ij}\) preserves the original structure of the equation, meaning that it involves the same number of layer potentials. For example, in the formulation using the weakly singular BM equation 3 , \(A_{ij}\) is expressed as \[\require{physics} \begin{align} A_{ij} = \\ \left( \begin{array}{c:c} D_{k_1}^{ij} + \delta_{ij} I/2 & -\varepsilon_1 S_{k_1}^{ij} \\ \hdashline&D_{k_0}^{ij} - \delta_{ij} I/2 \\ &+ \alpha \qty(N_{k_0}^{ij} - N_{k_1}^{ij})& \begin{align} &-\epsilon_0 S_{k_0}^{ij} \\ &- \alpha \qty(\epsilon_0 D_{k_0}^{*ij} - \epsilon_1 D_{k_1}^{*ij} ) \\ &- \alpha (\varepsilon_0 + \varepsilon_1) \delta_{ij} I/2 \end{align} \end{array} \right), \end{align}\] where \(Z_{k}^{ij}\) is the submatrix of a discretized layer potential \(Z_{k}\) associated with the index set \(I_i\) and \(I_j\), and \(\delta_{ij}\) is the Kronecker delta. Within each MPI process, the task of compressing the block linear system associated with the corresponding subtree is parallelized using OpenMP. By assigning the calculation of \(A_{ij}\) to a single OpenMP thread, the computational load can be evenly distributed. Note that \(A_{ij}\) is only calculated around the diagonal at the current level of the tree. To minimize MPI communication, each process locally stores one-dimensional data, such as global indices and the coordinates of the quadrature points. Within each assigned subtree, low-rank approximations are computed via the proxy-surface method. Because this method requires the skeleton indices [7] of adjacent cells in the entire tree, we perform a boundary exchange of these indices at the cells corresponding to the subtree endpoints, similar to the domain decomposition approach employed in finite difference methods.

4 Numerical experiments↩︎

We demonstrate the performance of the direct solver on the TSUBAME4.0 (T4) supercomputer at the Institute of Science Tokyo. As scattering obstacles, we employed the circular and star-shaped geometries illustrated in Fig. 3. For the circular geometry, the analytical solution to the transmission problem can be explicitly calculated. A plane wave propagating in the \(x_1\)-direction was used as the incident wave. The parameters used in the demonstration are listed in Table 1. The solver was primarily implemented in C++ using the Eigen (v3) library. In our parallel configuration, a single compute node (comprising 192 CPU cores) of T4 is assigned to each MPI process; thus, the number of nodes is identical to the number of processes.

We first verified the implementation for the circular scatterer using two compute nodes. As shown in Figure 4, we found that the direct solver based on the weakly singular BM equation 3 achieved \(O(h^3)\) convergence, while the solvers based on the other equations achieved only \(O(h)\). Here, \(h\) denotes the spacing between the quadrature nodes. In this figure and subsequent figures, we refer to the weakly singular BM equation as “WSBM.”

We next demonstrate the performance of the direct solver for each of the BIEs 1 –3 . Figure 5 illustrates the total computational time versus degrees of freedom using 1–8 compute nodes of T4. We observe that all the solvers exhibited nearly ideal weak and strong scalabilities. For example, for approximately \(1.31 \times 10^7\) degrees of freedom, the solver based on the weakly singular BM equation achieved a 7.7-fold speedup when using 8 nodes, corresponding to a strong scaling efficiency of 97%. The computational time for the standard BM equation was approximately 10%–20% shorter than that for the weakly singular BM equation. However, from the viewpoint of the convergence rate, the latter remains the most attractive formulation. We note that the computational time of the solver based on the PMCHWT equation 1 was almost identical to that for the weakly singular BM equation 3 , because the number and types of layer potentials are the same.

5 Conclusion↩︎

This paper described a scalable distributed parallel fast direct solver for the hierarchically semiseparable representation. The direct solver is based on the weakly singular BM equation with the one-point zeta-corrected quadrature and on the load-balancing technique accomplished by a block formulation of the linear system. Numerical results demonstrate that our direct solver achieves \(O(h^3)\) convergence and nearly ideal strong and weak scalabilities. As a future direction, it would be interesting to extend the solver to scale for even more practical problems by supporting a larger number of processes, including non-power-of-two configurations, and incorporating GPU acceleration.

Acknowledgment↩︎

This work was supported by the JSPS KAKENHI grant number 24K20783. It was carried out using the TSUBAME4.0 at the Institute of Science Tokyo, which was provided through the projects JHPCN and HPCI in Japan (project ID jh260048).

References↩︎