A Pre-splitting Green’s Function based Hybrid Fast Algorithmfor Multiscale Problems

Guang-Yu Zhu${}^{\mathbf{\text{1,2}}}$, Wei-Dong Li${}^{\mathbf{\text{2}}}$, Wei E. I. Sha${}^{\mathbf{\text{1}}}$, Hou-Xing Zhou${}^{\mathbf{\text{2}}}$, and Wei Hong${}^{\mathbf{\text{2}}}$

¹College of Information Science and Electronic Engineering
Zhejiang University, Hangzhou, 310027, China
gyzhu_china@163.com
²State Key Laboratory of Millimeter Waves
Frontiers Science Center for Mobile Information Communication and Security
Southeast University, Nanjing, 210096, China

Submitted On: February 19, 2023; Accepted On: July 26, 2023

ABSTRACT

Based on the splitting form of the Green’s function, a hybrid fast algorithm is proposed for efficient analysis of multiscale problems. In this algorithm, the Green’s function is a priori split into two parts: a spectrally band-limited part and a spatially localized part. Then, the fast Fourier transforms (FFT) utilizing the global Cartesian grid and the matrix compression method aided by an adaptive octree grouping are implemented for these two parts, respectively. Compared with the traditional methods which only employ the FFT for acceleration, the proposed hybrid fast algorithm is capable of maintaining low memory consumption in multiscale problems without compromising time cost. Moreover, the proposed algorithm does not need cumbersome geometric treatment to implement the hybridization, and can be established in a concise and straightforward manner. Several numerical examples discretized with multiscale meshes are provided to demonstrate the computational performance of proposed hybrid fast algorithm.

Index Terms: Fast algorithm, fast Fourier transform, Green’s function, matrix compression, multiscale problems.

I. INTRODUCTION

Many real-world electromagnetic (EM) problems require multiscale discretization [1, 2, 3, 4]. When the method of moments (MoM) [5], along with classical fast algorithms, is used to model the multi-scale problems, a single fast algorithm [6, 7, 8, 9, 10, 11, 12] is often insufficient to achieve satisfactory computational performance. Instead, due to capturing both the circuit physics and wave physics [1], hybrid fast algorithms [19, 20, 21, 22] generally provide a more practical path to efficiently solve EM multiscale problems.

As we know, the multilevel fast multipole algorithm (MLFMA) is widely used for efficient simulation of electrically large problems [6, 7, 8]. However, the MLFMA encounters the sub-wavelength breakdown [8] when dense mesh occurs. As a remedy, low-frequency fast algorithms [23, 24] were studied. Correspondingly, a series of hybrid algorithms were developed [25, 26, 27, 28, 29, 30], such as mixed-form FMA [19], interpolative decomposition (ID)-MLFMA [20], and adaptive cross approximation (ACA)-MLFMA [21].

Different from the MLFMA, the pre-correction-fast Fourier transforms (FFT)-based methods [9, 10, 11, 12, 13, 14, 15, 16, 17, 18] (e.g., adaptive integral method (AIM) [9], pre-corrected FFT (P-FFT) [10], and integral equation FFT (IE-FFT) [11]) are free from the sub-wavelength breakdown. Nevertheless, due to the intrinsic uniformity of the global Cartesian grids, the pre-correction-FFT-based methods do not allow for sufficient and varying spatial resolutions to handle multiscale discretizations. That is, by simply tuning the grid spacing, it is hard to achieve high computational efficiency for far interactions and maintain low storage efficiency for the near matrix at the same time [31].

Recently, in order to overcome this shortcoming of the pre-correction-FFT-based methods, a hybrid fast algorithm was developed in [31]. In this algorithm, the pre-correction-FFT-based methods is augmented with the matrix compression method (ACA) [36, 37, 38]. However, to implement such a hybridization, the procedures are not as straightforward as they seem [31]. In particular, a complicated spatial grouping strategy is inherently needed, and the relations between the required auxiliary geometric data structures (such as the Cartesian grid, the expansion box, and the tree-cube) are complex and require cautious treatments [31]. Such a situation motivates us to reconsider the overall framework of the pre-correction-FFT-based methods and wonder if there exists a more elegant way to develop a hybrid algorithm that consists of both the FFT and the matrix compression methods(e.g., ACA).

As an interesting alternative of the well-studied pre-correction-FFT-based methods, some pre-splitting-FFT-based method was developed in [32, 33]. Basically, the pre-splitting-FFT-based method bears similarities with the well-known pre-correction-FFT-based methods, since it also uses the Cartesian grid and the FFT as the essentials for the acceleration of the spatial convolutions. However, different from the widely-used pre-correction-FFT-based methods which typically utilize the Cartesian grid-based near-field corrections (based on some local expansion boxes) to compensate for the near-field contributions, the pre-splitting-FFT-based method [32, 33] instead resorts to some splitting form of the Green’s function so as to realize accurate and efficient computations. Despite of its unique feature, this pre-splitting-based framework has not attracted enough attentions in the computational electromagnetics (CEM) community, and a hybrid fast algorithm based on such framework has never been explored yet.

In this work, based on the splitting form of the Green’s function, a hybrid fast algorithm is proposed for efficiently analyzing multiscale problems. Specifically, the pre-splitting-FFT-based method originally devised for the problems with quasi-uniform discretizations is here further enhanced with the matrix compression method (ACA). Thus, a hybrid fast algorithm using both the FFT and the matrix compression method is established. In particular, the proposed hybrid fast algorithm is capable of maintaining low memory consumption in multiscale problems without compromising time cost, thus manifesting itself as a favorable multiscale extension of the pure FFT-based method. Furthermore, the required auxiliary geometric data structures herein (including the Cartesian grid and the spatial grouping octree) can be constructed independent of each other, and are not intertwined as in [31]. Consequently, the proposed hybrid algorithm can be implemented in a straightforward manner without any cumbersome geometric treatment.

In the following, after introducing the basic formulations and implementation key points in Section 2, several numerical examples discretized with multiscale meshes are then provided in Section 3 to demonstrate the computational performances of the proposed hybrid fast algorithm for multiscale problems.

II. FORMULATION

A perfectly electric conductor (PEC) object illuminated by an incident plane wave is considered. The combined field integral equation (CFIE) is employed to model this typical electromagnetic scattering problem. Throughout the paper, we use $\lambda$ to denote the wavelength in free space. Moreover, we use $k$ and $\eta$ to denote the wave number and the wave impedance in free space, respectively.

A. Combined field integral equation

The CFIE is the linear combination of the electric field integral equation (EFIE) and the magnetic field integral equation (MFIE), which is conventionally expressed as

$$\text{CFIE}=\alpha \text{EFIE}+\left(1-\alpha \right)\eta \text{MFIE},$$

(1)

where $\alpha$ denotes the combination factor and $0⩽\alpha ⩽1$, the EFIE is formulated as

$$\begin{array}{c}{{\mathbf{E}}^{\text{inc}}\left(\mathit{\text{r}}\right)|}_{\tan }=\hfill \\ \hfill {\left[jk\eta {\int }_S\left(\mathbf{J}\left({\mathit{\text{r}}}^{\prime }\right)+\frac{\nabla {\nabla }^{\prime }\cdot }{k^2}\mathbf{J}\left({\mathit{\text{r}}}^{\prime }\right)\right)\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\mathit{d}{s}^{\prime }\right]|}_{\tan },\end{array}$$

(2)

and the MFIE is formulated as

$$\begin{array}{c}{{\mathbf{H}}^{\text{inc}}\left(\mathit{\text{r}}\right)|}_{\tan }=\hfill \\ \hfill [\frac{\widehat{\mathit{\text{n}}}\times \mathbf{J}\left({\mathit{\text{r}}}^{\prime }\right)}{2}+P.V.{\int }_S\mathbf{J}\left({\mathit{\text{r}}}^{\prime }\right)\times \nabla \mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })d{s}^{\prime }]{|}_{\tan }.\end{array}$$

(3)

In the above, ${\mathbf{E}}^{\text{inc}}$ and ${\mathbf{H}}^{\text{inc}}$ are the incident electric field and the incident magnetic field, respectively. Moreover, $\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ denotes the Green’s function with the field point r and the source point ${\mathit{\text{r}}}^{\prime }$. Furthermore, $\mathbf{J}\left({\mathit{\text{r}}}^{\prime }\right)$ denotes the surface current.

By using the method of moments (MoM) [5] with the RWG function [39], the CFIE can be discretized into a matrix equation

$${\mathit{\text{Z}}}^{\left[N\times N\right]}\cdot {\mathit{\text{I}}}^{\left[N\times 1\right]}={\mathit{\text{V}}}^{\left[N\times 1\right]},$$

(4)

where Z is the impedance matrix, I is the unknown surface current coefficients vector, V is the righthand side vector of the incident field, and $N$ is the number of unknowns. More specifically, the impedance matrix Z corresponding to the CFIE is described by

$${\mathit{\text{Z}}}^{\text{CFIE}}=\alpha {\mathit{\text{Z}}}^{\text{EFIE}}+\left(1-\alpha \right)\eta {\mathit{\text{Z}}}^{\text{MFIE}},$$

(5)

where

$$\begin{array}{c}\hfill {\mathit{\text{Z}}}_{mn}^{\text{EFIE}}=jk\eta \underset{S_m}{\int }\mathit{d}s{\mathbf{f}}_m\left(\mathit{\text{r}}\right)\cdot \underset{S_n}{\int }\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime }){\mathbf{f}}_n\left({\mathit{\text{r}}}^{\prime }\right)\mathit{d}{s}^{\prime }\hfill \\ \hfill -j\frac{\eta }{k}\underset{S_m}{\int }\mathit{d}s\left[\nabla \cdot {\mathbf{f}}_m\left(\mathit{\text{r}}\right)\right]\cdot \underset{S_n}{\int }\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\left[\nabla \cdot {\mathbf{f}}_n\left({\mathit{\text{r}}}^{\prime }\right)\right]\mathit{d}{s}^{\prime },\hfill \end{array}$$

(6)

and

$$\begin{array}{c}\hfill {\mathit{\text{Z}}}_{mn}^{\text{MFIE}}=\frac{1}{2}\underset{S_m}{\int }{\mathbf{f}}_m\left(\mathit{\text{r}}\right)\cdot {\mathbf{f}}_n\left(\mathit{\text{r}}\right)\mathit{d}s\hfill \\ \hfill +\underset{S_m}{\int }\mathit{d}s\left[\widehat{\mathit{\text{n}}}\times {\mathbf{f}}_m\left(\mathit{\text{r}}\right)\right]\cdot \underset{S_n}{\int }\nabla \mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\times {\mathbf{f}}_n\left({\mathit{\text{r}}}^{\prime }\right)\mathit{d}{s}^{\prime }.\hfill \end{array}$$

(7)

In the above, $m=1,\mathrm{\cdots },N$ and $n=1,\mathrm{\cdots },N$. Moreover, ${\mathbf{f}}_m$ and ${\mathbf{f}}_n$ denote the $m$th and $n$th RWG functions, respectively. Furthermore, $S_m$ and $S_n$ are the supports of the corresponding RWG functions, respectively.

B. Pre-splitting the Green’s function

In the following, the splitting forms of the Green’s function and its gradient are introduced. First, consider the Green’s function in the EFIE, namely

$$\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })=\frac{e^{-jkR}}{4\pi R},$$

(8)

where

$$\begin{array}{c}\hfill R=\left|\mathit{\text{r}}-{\mathit{\text{r}}}^{\prime }\right|\hfill \\ \hfill =\sqrt{{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2+{\left(z-z^{\prime }\right)}^2}\hfill \end{array}$$

(9)

denotes the distance between the field point $\mathbf{r}$ and the source point ${\mathbf{r}}^{\prime }$. Clearly, the above can also be equivalently represented as

$$\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })=\mathrm{Re}\left[\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]+j\mathrm{Im}\left[\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right],$$

(10)

where

$$\mathrm{Re}\left[\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]=\frac{\cos kR}{4\pi R},$$

(11)

$$\mathrm{Im}\left[\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]=-\frac{\sin kR}{4\pi R}.$$

(12)

A basic property of $\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ is that, this function is singular for $R=0$, and its magnitude is globally nonzero for $R>0$. Also, note that the singular behavior of $\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ comes from $\mathrm{Re}\left[\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]$.

Using the pre-splitting strategy [32, 33], the Green’s function (8) can be written as two parts, namely

$$\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })={\mathrm{G}}_{\text{E}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })+{\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime }),$$

(13)

where

$${\mathrm{G}}_{\text{E}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })=\{\begin{array}{cccccccccccccccccccc}\hfill Re\left[\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]-\mathrm{\Phi }\left(R\right),R<\delta \hfill & & & & & & & & & & & & & & & & & & & \\ \hfill 0\begin{array}{}\end{array},R\ge \delta \hfill & & & & & & & & & & & & & & & & & & & \end{array},$$

(14)

$${\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })=\{\begin{array}{cccccccccccccccccccc}\hfill \mathrm{\Phi }\left(R\right)+jIm\left[\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right],R<\delta \hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\begin{array}{}\end{array},R\ge \delta \hfill & & & & & & & & & & & & & & & & & & & \end{array},$$

(15)

with $\delta$ being a splitting threshold, and $\mathrm{\Phi }\left(R\right)$ being some auxiliary function. Particularly, the above two parts manifest the following properties.

• Function ${\mathrm{G}}_{\text{E}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ is a singular and spatially localized function.

It has nonzero function value only when $R<\delta$, and ${\mathrm{G}}_E(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })=0$ for $R\ge \delta$.

• Function ${\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ is a globally oscillatory and spectrally band-limited function.

It is smooth and bounded for $\forall \mathit{\text{r}}$ and $\forall {\mathit{\text{r}}}^{\prime }$.

To have more intuitive understanding on the properties of function ${\mathrm{G}}_{\text{E}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ and ${\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$, the corresponding function images are shown in Fig. 1.

Moreover, from a physical and spectral perspective, ${\mathrm{G}}_{\text{E}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ mainly captures the evanescent waves, and ${\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ mainly captures the propagating waves.

Figure 1: The Green’s function is a priori split into a singular but spatially localized part ${\mathrm{G}}_E$ and a globally oscillatory but spectrally band-limited part ${\mathrm{G}}_P$. The relevant magnitudes of the original Green’s function and of these two parts are visualized in the figure.

To achieve the above splitting and corresponding properties, the auxiliary function $\mathrm{\Phi }\left(R\right)$ is specially designed and one simple way is to choose $\mathrm{\Phi }\left(R\right)$ to be a polynomial of the following form [32, 33]

$$\mathrm{\Phi }\left(R\right)=a{R}^3+b{R}^2+c,$$

(16)

where $a$, $b$, and $c$ are unknown coefficients. Clearly, the above is a three order polynomial without the linear term, and it implies

$$d\mathrm{\Phi }\left(0\right)/dR=0,$$

(17)

meaning that the first order derivative function of $\mathrm{\Phi }\left(R\right)$ is zero at $R=0$ (where $\mathit{\text{r}}={\mathit{\text{r}}}^{\prime }$).

Given the splitting threshold $\delta$, the coefficients in (16) can be determined by some continuity conditions at $R=\delta$ as follows

$$\{\begin{array}{cccccccccccccccccccc}\mathrm{\Phi }\left(\delta \right)\hfill & =\hfill & \mathrm{Re}\left[\mathrm{G}\left(\delta \right)\right]\hfill & & & & & & & & & & & & & & & & & \\ d\mathrm{\Phi }\left(\delta \right)/dR\hfill & =\hfill & d\mathrm{Re}\left[\mathrm{G}\left(\delta \right)\right]/dR\hfill & & & & & & & & & & & & & & & & & \\ d^2\mathrm{\Phi }\left(\delta \right)/d{R}^2\hfill & =\hfill & d^2\mathrm{Re}\left[\mathrm{G}\left(\delta \right)\right]/d{R}^2\hfill & & & & & & & & & & & & & & & & & \end{array},$$

(18)

or, by solving the following matrix equation

$$\left[\begin{array}{ccc}\hfill {\delta }^3\hfill & \hfill {\delta }^2\hfill & \hfill 1\hfill \\ \hfill 3{\delta }^2\hfill & \hfill 2\delta \hfill & \hfill 0\hfill \\ \hfill 6\delta \hfill & \hfill 2\hfill & \hfill 0\hfill \end{array}\right]\left[\begin{array}{cccccccccccccccccccc}\hfill a\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill b\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill c\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right]=\left[\begin{array}{cccccccccccccccccccc}\mathrm{Re}\left[\mathrm{G}\left(\delta \right)\right]\hfill & & & & & & & & & & & & & & & & & & & \\ d\mathrm{Re}\left[\mathrm{G}\left(\delta \right)\right]/dR\hfill & & & & & & & & & & & & & & & & & & & \\ d^2\mathrm{Re}\left[\mathrm{G}\left(\delta \right)\right]/d{R}^2\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right].$$

(19)

The derivative functions of $\mathrm{Re}\left[\mathrm{G}\left(R\right)\right]$ used in the above are given in the appendix. Thus, after $\mathrm{\Phi }\left(R\right)$ is obtained, the splitting representation (13) is completely determined.

The polynomial choice (16), the implicit continuity condition (17), and the explicit continuity condition (18) are based on the consideration that, after the splitting, the resultant ${\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ (or, more precisely, $\mathrm{Re}\left[{\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]$) should be smooth everywhere for $\forall \mathit{\text{r}}$ and $\forall {\mathit{\text{r}}}^{\prime }$, so that it can be interpolated accurately without much difficulty (the interpolation procedures will be discussed in the next subsection).

To make it more clear, without loss of generality, consider the following representative situation. The source point ${\mathit{\text{r}}}^{\prime }$ is fixed to the coordinate origin, and the field point r is moving along the $x$ axis. Then, the corresponding function values of $\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ and ${\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ are evaluated, as shown in Fig. 3. We can clearly see that $\mathrm{Re}\left[\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]$ is singular when r approaches ${\mathit{\text{r}}}^{\prime }$, while $\mathrm{Re}\left[{\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]$ is continuous and smooth everywhere. Here, $\mathrm{Re}\left[{\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]$ is smooth at $R=0$ (where $\mathit{\text{r}}={\mathit{\text{r}}}^{\prime }$) as a result of (17), and smooth at the splitting point $R=\delta$ as a result of (18).

Next, consider the gradient of the Green’s function in the MFIE, namely

$$\nabla \mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })=\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(R\right)\nabla R,$$

(20)

with

$$\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(R\right)=\mathrm{Re}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(R\right)\right]+j\mathrm{Im}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(R\right)\right],$$

(21)

where

$$\mathrm{Re}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(R\right)\right]=-\frac{\cos kR+kR\sin kR}{4\pi R^2},$$

(22)

$$\mathrm{Im}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(R\right)\right]=-\frac{kR\cos kR-\sin kR}{4\pi R^2},$$

(23)

and

$$\nabla R=\frac{x-x^{\prime }}{R}\widehat{\mathit{\text{x}}}+\frac{y-y^{\prime }}{R}\widehat{\mathit{\text{y}}}+\frac{z-z^{\prime }}{R}\widehat{\mathit{\text{z}}}.$$

(24)

Similar to (8), function $\nabla \mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ is singular for $R=0$, and its magnitude is globally nonzero for $R>0$. Also, note that the singular behavior of $\nabla \mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ comes from $\mathrm{Re}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]$.

Again, using the pre-splitting strategy [32, 33], the gradient Green’s function (20) can be written as two parts, namely

$$\nabla \mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })=\nabla {\mathrm{G}}_{\text{E}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })+\nabla {\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime }),$$

(25)

where

$$\nabla {\mathrm{G}}_{\text{E}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })=\{\begin{array}{cccccccccccccccccccc}\hfill \left\{\mathrm{Re}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(R\right)\right]-{\mathrm{\Phi }}_{\text{g}}\left(R\right)\right\}\nabla R,R<\delta \hfill & & & & & & & & & & & & & & & & & & & \\ \hfill 0\begin{array}{}\end{array},R⩾\delta \hfill & & & & & & & & & & & & & & & & & & & \end{array},$$

(26)

$$\nabla {\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })=\{\begin{array}{cccccccccccccccccccc}\hfill \left\{{\mathrm{\Phi }}_{\text{g}}\left(R\right)+j\text{Im}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(R\right)\right]\right\}\nabla R,R<\delta \hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \nabla \mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\begin{array}{}\end{array},R⩾\delta \hfill & & & & & & & & & & & & & & & & & & & \end{array},$$

(27)

with $\delta$ being the splitting threshold, and ${\mathrm{\Phi }}_{\text{g}}\left(R\right)$ being some auxiliary function.

Here, the following auxiliary function is used to achieve the desired splitting, namely

$${\mathrm{\Phi }}_{\text{g}}\left(R\right)=a{R}^3+b{R}^2+cR.$$

(28)

Note that, different from (16), the above is a three order polynomial without the constant term, and it implies

$${\mathrm{\Phi }}_{\text{g}}\left(0\right)=0,$$

(29)

which clearly means that the function value of ${\mathrm{\Phi }}_{\text{g}}\left(R\right)$ is zero at $R=0$ (where $\mathit{\text{r}}={\mathit{\text{r}}}^{\prime }$).

Given the splitting threshold $\delta$, the coefficients in (28) can be determined by some continuity conditions at $R=\delta$ as follows

$$\{\begin{array}{cccccccccccccccccccc}{\mathrm{\Phi }}_{\text{g}}\left(\delta \right)\hfill & =\hfill & \text{Re}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(\delta \right)\right]\hfill & & & & & & & & & & & & & & & & & \\ d{\mathrm{\Phi }}_{\text{g}}\left(\delta \right)/dR\hfill & =\hfill & d\text{Re}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(\delta \right)\right]/dR\hfill & & & & & & & & & & & & & & & & & \\ d^2{\mathrm{\Phi }}_{\text{g}}\left(\delta \right)/d{R}^2\hfill & =\hfill & d^2\text{Re}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(\delta \right)\right]/d{R}^2\hfill & & & & & & & & & & & & & & & & & \end{array},$$

(30)

or, by solving the following matrix equation

$$\left[\begin{array}{ccc}\hfill {\delta }^3\hfill & \hfill {\delta }^2\hfill & \hfill \delta \hfill \\ \hfill 3{\delta }^2\hfill & \hfill 2\delta \hfill & \hfill 1\hfill \\ \hfill 6\delta \hfill & \hfill 2\hfill & \hfill 0\hfill \end{array}\right]\left[\begin{array}{cccccccccccccccccccc}\hfill a\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill b\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill c\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right]=\left[\begin{array}{cccccccccccccccccccc}\mathrm{Re}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(\delta \right)\right]\hfill & & & & & & & & & & & & & & & & & & & \\ d\mathrm{Re}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(\delta \right)\right]/dR\hfill & & & & & & & & & & & & & & & & & & & \\ d^2\mathrm{Re}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(\delta \right)\right]/d{R}^2\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right].$$

(31)

The derivative functions of $\mathrm{Re}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(R\right)\right]$ used in the above are given in the appendix. Thus, after ${\mathrm{\Phi }}_{\text{g}}\left(R\right)$ is obtained, the splitting representation (25) is completely determined.

Similar as before, (28), (29), and (30) are based on the consideration that, after the splitting, the resultant $\nabla {\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ (or, more precisely, each Cartesian component $\mathrm{Re}{[\nabla {\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })]}_{\sigma }$, $\sigma =\{x,y,z\}$) should be smooth everywhere for $\forall \mathit{\text{r}}$ and $\forall {\mathit{\text{r}}}^{\prime }$, so that it can be interpolated accurately without much difficulty.

To make it more clear, consider the representative situation of ${\mathit{\text{r}}}^{\prime }$ and r as before. Then, the corresponding function values of ${\left[\nabla \mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]}_x$ and ${\left[\nabla {\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]}_x$ are evaluated, as shown in Fig. 4. Here, $\mathrm{Re}{[\nabla {\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })]}_x$ is smooth at $R=0$ (where $\mathit{\text{r}}={\mathit{\text{r}}}^{\prime }$) as a result of (29), and smooth at the splitting point $R=\delta$ as a result of (30).

Note that the auxiliary function (28) chosen for the splitting of $\mathrm{Re}{[\nabla \mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })]}_{\sigma }$ ($\sigma =\{x,y,z\}$) is different from the auxiliary function (16) chosen for the splitting of $\mathrm{Re}\left[\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]$. Such difference is due to the fact that, unlike $\mathrm{Re}\left[\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]$, $\mathrm{Re}{[\nabla \mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })]}_{\sigma }$ ($\sigma =\{x,y,z\}$) manifest different (i.e., positive or negative) singular behaviors when r approaches ${\mathit{\text{r}}}^{\prime }$ along different directions (e.g., along $+x$ or $-x$ direction), which can be deduced from the definition of (24), or, more intuitively, observed from Fig. 3 and Fig. 4.

Figure 2: Systematic diagram of the flow chart and relevant key elements of the presented hybrid fast algorithm. Particularly, the two auxiliary geometric data structures, including the Cartesian grid and the octree grouping, can be constructed independently.

Figure 3: The function values of $\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ and ${\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$. The source point ${\mathit{\text{r}}}^{\prime }$ is fixed to the coordinate origin, and the field point r is moving along the $x$ axis.

Figure 4: The function values of ${\left[\nabla \mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]}_x$ and ${\left[\nabla {\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })\right]}_x$. The source point ${\mathit{\text{r}}}^{\prime }$ is fixed to the coordinate origin, and the field point r is moving along the $x$ axis.

With the splitting form of the Green’s function (13) and its gradient (25), the impedance matrix ${\mathit{\text{Z}}}^{\text{CFIE}}$ can be correspondingly split into two parts, yielding

$${\mathit{\text{Z}}}^{\text{CFIE}}={\mathit{\text{Z}}}_{\text{E}}+{\mathit{\text{Z}}}_{\text{P}},$$

(32)

where ${\mathit{\text{Z}}}_{\text{E}}$ is related to ${\mathrm{G}}_{\text{E}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ and $\nabla {\mathrm{G}}_{\text{E}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$, while ${\mathit{\text{Z}}}_{\text{P}}$ is related to ${\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ and $\nabla {\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$. For convenience, ${\mathit{\text{Z}}}_{\text{E}}$ and ${\mathit{\text{Z}}}_{\text{P}}$ are termed here as the evanescent matrix and the propagating matrix, respectively. When the commonly-used Krylov subspace iterative method [35] is employed to solve the matrix equation (4), the matrix-vector product ${\mathit{\text{Z}}}^{\text{CFIE}}\cdot \mathit{\text{I}}$ is computed. Using the splitting representation (32), this product can be equivalently implemented as

$${\mathit{\text{Z}}}^{\text{CFIE}}\cdot \mathit{\text{I}}={\mathit{\text{Z}}}_{\text{E}}\cdot \mathit{\text{I}}+{\mathit{\text{Z}}}_{\text{P}}\cdot \mathit{\text{I}}.$$

(33)

As will be shown in the following subsections, by fully exploiting the properties of the splitting form of the Green’s function and its gradient, both ${\mathit{\text{Z}}}_{\text{E}}\cdot \mathit{\text{I}}$ and ${\mathit{\text{Z}}}_{\text{P}}\cdot \mathit{\text{I}}$ can be evaluated in an efficient manner. For clarity, before going into further discussion, basic concepts and key elements of the presented hybrid fast algorithm are illustrated in Fig. 2.

C. Propagating matrix calculation

Due to the unique properties (smoothness and translation invariance) of ${\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ and $\nabla {\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$, the matrix-vector product ${\mathit{\text{Z}}}_{\text{P}}\cdot \mathit{\text{I}}$ can be calculated accurately and efficiently, with the aid of proper interpolation and theFFT.

To this end, some auxiliary geometric data structures need to be defined first. Specifically, the object considered is enclosed by a three dimensional (3D) rectangular bounding box. Then, this bounding box is subdivided regularly along each Cartesian dimension, yielding a cluster of global Cartesian grids [9, 10, 11, 12, 13, 14]. The grid spacings along these Cartesian dimensions are here denoted by $h_x$, $h_y$, $h_z$, respectively. The number of grids along these Cartesian dimensions are denoted by $N_x$, $N_x$, $N_x$, respectively. Thus, a total of $N_{\text{Grid}}=N_x\times N_y\times N_z$ Cartesian grids are defined.

Recall that ${\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ and $\nabla {\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ obtained from the splitting process are smooth everywhere for $\forall \mathit{\text{r}}$ and $\forall {\mathit{\text{r}}}^{\prime }$. Thus, they can be well approximated using the following Cartesian grid-based sampling expansions

$${\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })=\sum _{\mathit{\text{u}}\in B_m}\sum _{\mathit{\text{v}}\in B_n}{\beta }_{\mathit{\text{u}}}\left(\mathit{\text{r}}\right){\mathrm{G}}_{\text{P}}(\mathit{\text{u}},\mathit{\text{v}}){\beta }_{\mathit{\text{v}}}\left({\mathit{\text{r}}}^{\prime }\right),$$

(34)

$$\nabla {\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })=\sum _{\mathbf{u}\in B_m}\sum _{\mathit{\text{v}}\in B_n}{\beta }_{\mathit{\text{u}}}\left(\mathit{\text{r}}\right)\nabla {\mathrm{G}}_{\text{P}}(\mathit{\text{u}},\mathit{\text{v}}){\beta }_{\mathit{\text{v}}}\left({\mathit{\text{r}}}^{\prime }\right),$$

(35)

for $\forall \mathit{\text{r}}$ and $\forall {\mathit{\text{r}}}^{\prime }$, where $B_m$ and $B_n$ denote the corresponding local expansion boxes that enclose r and ${\mathit{\text{r}}}^{\prime }$, respectively. Note that, depending on the relative position between r and ${\mathit{\text{r}}}^{\prime }$, $B_m$ and $B_n$ can be either well-separated or overlapping, as illustrated in Fig. 5. Moreover, u and v denote the coordinates of the Cartesian grids related to $B_m$ and $B_n$, respectively. Furthermore, ${\beta }_{\mathit{\text{u}}}\left(\mathit{\text{r}}\right)$ (or ${\beta }_{\mathit{\text{v}}}\left({\mathit{\text{r}}}^{\prime }\right)$) denotes the 3D Lagrange interpolation basis function based on the Cartesian grids of $B_m$ (or $B_n$), which is defined concretely by

$${\beta }_{\mathit{\text{u}}}\left(\mathit{\text{r}}\right)={\beta }_{p_x}\left(x\right){\beta }_{p_y}\left(y\right){\beta }_{p_z}\left(z\right),$$

(36)

where

$${\beta }_{p_{\sigma }}\left(\sigma \right)=\prod _{p_{\sigma }=0,p_{\sigma }\ne p^{\prime }{}_{\sigma }}^{M_{\sigma }}\frac{\left(\sigma -{\sigma }_{p^{\prime }{}_{\sigma }}\right)}{\left({\sigma }_{p_{\sigma }}-{\sigma }_{p^{\prime }{}_{\sigma }}\right)},\sigma =\{x,y,z\},$$

(37)

with $M_{\sigma }$ being the expansion order and ${\sigma }_{p_{\sigma }}$ being the coordinates of the Cartesian grids involved. Clearly, when $M_x=M_y=M_z=M$, a total of ${\left(M+1\right)}^3$ Cartesian grids are involved in (36).

Figure 5: Computation corresponding to the propagating matrix ${\mathit{\text{Z}}}_{\text{P}}$ can be accelerated by the FFT using the Global Cartesian grid and the local expansion boxes. Typical examples of the expansion boxes are illustrated using different colors. Here, Box 1 and Box 2 are well-separated, while Box 2 and Box 3 are overlapping.

With the expansions (34) and (35), ${\mathbf{Z}}_{\text{P}}^{\text{CFIE}}\cdot \mathit{\text{I}}$ can be correspondingly rewritten as

$$\begin{array}{c}\hfill {\mathit{\text{Z}}}_{\text{P}}^{\text{CFIE}}\cdot \mathit{\text{I}}=j\mathrm{\Phi }k\eta \overline{\mathrm{\Pi }}\cdot \left[{\mathrm{G}}_{\text{P}}(\mathit{\text{u}},\mathit{\text{v}})\right]\cdot {\overline{\mathrm{\Pi }}}^T\cdot \mathit{\text{I}}\hfill \\ \hfill -j\alpha \frac{\eta }{k}{\mathrm{\Pi }}_{\text{d}}\cdot \left[{\mathrm{G}}_{\text{P}}(\mathit{\text{u}},\mathit{\text{v}})\right]\cdot {\mathrm{\Pi }}_d^T\cdot \mathit{\text{I}}\hfill \\ \hfill +\left(1-\alpha \right)\eta {\overline{\mathrm{\Pi }}}_{\text{g}}\cdot \left[\nabla {\mathrm{G}}_{\text{P}}(\mathit{\text{u}},\mathit{\text{v}})\right]\times {\overline{\mathrm{\Pi }}}^T\cdot \mathit{\text{I}},\hfill \end{array}$$

(38)

where $\left[{\mathrm{G}}_{\text{P}}\right]$ and ${\left[\nabla {\mathrm{G}}_{\text{P}}\right]}_{\sigma }$ are triple block Toeplitz matrices of dimension $N_{\text{Grid}}\times N_{\text{Grid}}$. Moreover, $\overline{\mathrm{\Pi }}$, ${\mathrm{\Pi }}_{\text{d}}$, and ${\overline{\mathrm{\Pi }}}_{\text{g}}$ are sparse matrices of dimension $N\times N_{\text{Grid}}$, defined concretely by

$$\overline{\mathrm{\Pi }}=\underset{S}{\int }\left[\begin{array}{cccccccccccccccccccc}\hfill {\mathbf{f}}_1\left(\mathit{\text{r}}\right)\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mathbf{f}}_2\left(\mathit{\text{r}}\right)\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mathbf{f}}_N\left(\mathit{\text{r}}\right)\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right]\left[\begin{array}{cccccccccccccccccccc}\hfill {\beta }_1\left(\mathit{\text{r}}\right),\hfill & \hfill \mathrm{\cdots }\hfill & \hfill ,{\beta }_{N_{\text{Grid}}}\left(\mathit{\text{r}}\right)\hfill & & & & & & & & & & & & & & & & & \end{array}\right]\mathit{d}s,$$

(39)

$${\mathrm{\Pi }}_{\text{d}}=\underset{S}{\int }\left[\begin{array}{cccccccccccccccccccc}\hfill \nabla \cdot {\mathbf{f}}_1\left(\mathit{\text{r}}\right)\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \nabla \cdot {\mathbf{f}}_2\left(\mathit{\text{r}}\right)\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \nabla \cdot {\mathbf{f}}_N\left(\mathit{\text{r}}\right)\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right]\left[\begin{array}{cccccccccccccccccccc}\hfill {\beta }_1\left(\mathit{\text{r}}\right),\hfill & \hfill \mathrm{\cdots }\hfill & \hfill ,{\beta }_{N_{\text{Grid}}}\left(\mathit{\text{r}}\right)\hfill & & & & & & & & & & & & & & & & & \end{array}\right]\mathit{d}s,$$

(40)

$${\overline{\mathrm{\Pi }}}_{\text{g}}=\underset{S}{\int }\left[\begin{array}{cccccccccccccccccccc}\hfill \widehat{\mathbf{n}}\times {\mathbf{f}}_1\left(\mathit{\text{r}}\right)\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \widehat{\mathbf{n}}\times {\mathbf{f}}_2\left(\mathit{\text{r}}\right)\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \widehat{\mathbf{n}}\times {\mathbf{f}}_N\left(\mathit{\text{r}}\right)\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right]\left[\begin{array}{cccccccccccccccccccc}\hfill {\beta }_1\left(\mathit{\text{r}}\right),\hfill & \hfill \mathrm{\cdots }\hfill & \hfill ,{\beta }_{N_{\text{Grid}}}\left(\mathit{\text{r}}\right)\hfill & & & & & & & & & & & & & & & & & \end{array}\right]\mathit{d}s.$$

(41)

Exploiting the block Toeplitz structure of $\left[{\mathrm{G}}_{\text{P}}\right]$ and ${\left[\nabla {\mathrm{G}}_{\text{P}}\right]}_{\sigma }$, the matrix-vector product ${\mathit{\text{Z}}}_P\cdot \mathit{\text{I}}$ can then be evaluated efficiently by the FFT as follows

$$\begin{array}{c}\hfill {\mathit{\text{Z}}}_{\text{P}}^{\text{CFIE}}\cdot \mathit{\text{I}}=j\alpha k\eta \overline{\mathrm{\Pi }}\cdot {ℱ}^{-1}\left\{ℱ\left\{\left[{\mathrm{G}}_{\text{P}}\right]\right\}\cdot ℱ\left\{{\overline{\mathrm{\Pi }}}^T\cdot \mathbf{I}\right\}\right\}\hfill \\ \hfill -j\alpha \frac{\eta }{k}{\overline{\mathrm{\Pi }}}_{\text{d}}\cdot {ℱ}^{-1}\left\{ℱ\left\{\left[{\mathrm{G}}_{\text{P}}\right]\right\}\cdot ℱ\left\{{\overline{\mathrm{\Pi }}}_d^T\cdot \mathbf{I}\right\}\right\}\hfill \\ \hfill +\left(1-\alpha \right)\eta {\overline{\mathrm{\Pi }}}_{\text{g}}\cdot {ℱ}^{-1}\left\{ℱ\left\{\left[\nabla {\mathrm{G}}_{\text{P}}\right]\right\}\times ℱ\left\{{\overline{\mathrm{\Pi }}}^T\cdot \mathbf{I}\right\}\right\},\hfill \end{array}$$

(42)

where $ℱ$ and ${ℱ}^{-1}$ denote the forward and inverse FFT.

Notably, in the conventional pre-correction-based framework [9, 10, 11, 12, 13], the original Green’s function $\mathrm{G}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ is approximated with the grid-based sampling expansion. Differently, for the pre-splitting-based framework [32, 33] adopted herein, the function ${\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ is instead approximated with the grid-based sampling expansion, as in (34). Since ${\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ is smooth everywhere without singularity, it is possible to approximate ${\mathrm{G}}_{\text{P}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ with sufficiently good accuracy for arbitrary r and ${\mathit{\text{r}}}^{\prime }$ by proper interpolation. Consequently, unlike [9, 10, 11, 12, 13, 31], the Cartesian grid-based near-field corrections are not necessary here.

D. Evanescent matrix calculation

Referring to (14) and (26), we know that ${\mathit{\text{Z}}}_{\text{E}}$, dominated by ${\mathrm{G}}_{\text{E}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$ and $\nabla {\mathrm{G}}_{\text{E}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$, has nonzero elements only when $R<\delta$. Therefore, with quasi-uniform meshes, ${\mathit{\text{Z}}}_{\text{E}}$ is a typical sparse matrix with few nonzero elements, and thus it can be stored without much effort. However, in multiscale problems, dense mesh always occurs. Thus, near-field interactions within $R<\delta$ increase substantially. In this case, the nonzero elements of ${\mathit{\text{Z}}}_{\text{E}}$ will increase dramatically.

Fortunately, in the pre-splitting-based framework adopted, this difficulty can be trivially addressed. In a nut shell, we first construct a spatial grouping only based on the geometric position of the basis functions, and then apply the matrix compression methods based on such a grouping to achieve memory reductions of ${\mathit{\text{Z}}}_{\text{E}}$.

Here, a standard adaptive octree is utilized for accomplishing the spatial grouping, and the basic procedures for its construction are outlined as follows.

1. A root cube enclosing the whole object is constructed.

2. The nonempty cubes are subdivided into eight sub-cubes recursively until the number of basis functions within each cube is smaller than a prescribed constant.

3. An octree is built when the process finishes.

Accordingly, several concepts related to the octree are then defined as follows.

•Those cubes that do not need to be further subdivided are called leaf cubes. Here, the leaf cubes are denoted by $C_s$, with $s$ being the index of some leaf cube. It should be noted that not all of the leaf cubes are located at the same octree level. Typical examples of the leaf cubes are shown in Fig. 6 and are highlighted with colors.

• Two leaf cubes (e.g., $C_s$ and $C_t$) are called neighbor cubes if they share at least one vertex, such as Cube 1 and Cube 5; otherwise, they are called well-separated cubes, such as Cube 1 and Cube 3. For convenience, we use $C_s\cap C_t\ne \mathrm{\varnothing }$ to denote that two cubes are neighbor cubes, and use $C_s\cap C_t=\mathrm{\varnothing }$ to denote that two cubes are well-separated cubes.

Note that, in [31], the conventional octree grouping cannot be directly applied, and a more involved grouping scheme is developed, in which the cubes for grouping the basis functions are induced from the Cartesian grids. In contrast, here the basis function grouping does not have to be dependent on the Cartesian grids, and can be implemented independently and conventionally. In other words, the commonly-used adaptive octree can be directly employed to achieve the required groupings, and the treatments of the required auxiliary geometric data structures herein are thus simple and straightforward.

Figure 6: The storage of ${\mathit{\text{Z}}}_E$ is effectively compressed by using the spatial grouping octree and the matrix compression method. Illustration of the octree grouping is shown in the left subfigure, where some leaf cubes are highlighted using colors. Here, Cube 5 is neighbor cube of Cube 1, while Cube 2, Cube 3, Cube 4 are all well-separated cubes of Cube 1. Besides, the nonzero range of ${\mathrm{G}}_{\text{E}}(\mathit{\text{r}},{\mathit{\text{r}}}^{\prime })$, characterized by $R<\delta$, is illustrated in the right subfigure. Note that the establishment of the grouping and the application of the compression can be realized independent of the Cartesian grids and the FFT.

On the basis of the spatial groupings above, ${\mathit{\text{Z}}}_{\text{E}}$ can then be rewritten as

$${\mathit{\text{Z}}}_{\text{E}}={\left\{{\mathit{\text{Z}}}_{\text{E}}^{\mathit{\text{st}}}\right\}}_{C_s\cap C_t\ne \mathrm{\varnothing }}+{\left\{{\mathit{\text{Z}}}_{\text{E}}^{\mathit{\text{st}}}\right\}}_{C_s\cap C_t=\mathrm{\varnothing }},$$

(43)

where

• ${\left\{{\mathit{\text{Z}}}_{\text{E}}^{\mathit{\text{st}}}\right\}}_{C_s\cap C_t\ne \mathrm{\varnothing }}$ denote those matrix subblocks corresponding to the interactions between neighbor cubes. They are directly computed and stored.

• ${\left\{{\mathit{\text{Z}}}_{\text{E}}^{\mathit{\text{st}}}\right\}}_{C_s\cap C_t=\mathrm{\varnothing }}$ denote those matrix subblocks corresponding to the interactions between well-separated cubes. These matrix subblocks can be low-rank. They can be compressed and stored in a compact form.

Here, the ACA [36, 37, 38] is employed for accomplishing the matrix compression. Detailed procedures of the ACA can be found in, for example, [37, 38]. Correspondingly, ${\left\{{\mathit{\text{Z}}}_{\text{E}}^{\mathit{\text{st}}}\right\}}_{C_s\cap C_t=\mathrm{\varnothing }}$ can be then factorized as

$${\left\{{\mathit{\text{Z}}}_{\text{E}}^{\mathit{\text{st}}}\right\}}_{C_s\cap C_t=\mathrm{\varnothing }}={\left\{{\mathit{\text{U}}}_{\text{E}}^{\mathit{\text{s}}}\cdot {\mathit{\text{V}}}_{\text{E}}^{\mathit{\text{t}}}\right\}}_{C_s\cap C_t=\mathrm{\varnothing }}.$$

(44)

Thus, the matrix-vector product ${\mathit{\text{Z}}}_{\text{E}}\cdot \mathit{\text{I}}$ can be calculated as

$$\begin{array}{c}\hfill {\mathit{\text{Z}}}_{\text{E}}\cdot \mathit{\text{I}}={\left\{{\mathit{\text{Z}}}_{\text{E}}^{\mathit{\text{st}}}\right\}}_{C_s\cap C_t\ne \mathrm{\varnothing }}\cdot \mathit{\text{I}}+{\left\{{\mathit{\text{Z}}}_{\text{E}}^{\mathbf{st}}\right\}}_{C_s\cap C_t=\mathrm{\varnothing }}\cdot \mathit{\text{I}}\hfill \\ \hfill ={\left\{{\mathit{\text{Z}}}_{\text{E}}^{\mathit{\text{st}}}\cdot {\mathit{\text{I}}}^t\right\}}_{C_s\cap C_t\ne \mathrm{\varnothing }}+{\left\{{\mathit{\text{U}}}_{\text{E}}^{\mathit{\text{s}}}\cdot \left({\mathit{\text{V}}}_{\text{E}}^{\mathbf{t}}\cdot {\mathit{\text{I}}}^t\right)\right\}}_{C_s\cap C_t=\mathrm{\varnothing }}.\hfill \end{array}$$

(45)

Some implementation details for the compression of ${\left\{{\mathbf{Z}}_{\text{E}}^{\mathit{\text{st}}}\right\}}_{C_s\cap C_t=\mathrm{\varnothing }}$ are further clarified as follows. In particular, as shown in Fig. 6, due to the localized nature of ${\mathit{\text{Z}}}_{\text{E}}$ decided by $R<\delta$, for a basis function $\varphi$ in Cube 1, all the basis functions in Cube 2 have full interactions with $\varphi$, while all the basis functions in Cube 4 have zero interactions with $\varphi$. Meanwhile, only parts of the basis functions in Cube 3 and Cube 5 have nonzero interactions with $\varphi$. Hence, for two well separated groups with index sets s and t, only the submatrices representing nonzero interactions is to be compressed. Thus, the matrix ${\mathit{\text{Z}}}_{\text{E}}^{\mathit{\text{st}}}$ which represents the interactions between the two groups is in fact described more precisely by

$${\mathit{\text{Z}}}_{\text{E}}^{\mathit{\text{st}}}={\mathit{\text{P}}}^{\mathit{\text{s}}\mathit{\text{s}}}\left[\begin{array}{cccccccccccccccccccc}\hfill {\mathit{\text{Z}}}_{\text{E}}^{{\mathit{\text{s}}}^{\prime }{\mathit{\text{t}}}^{\prime }}\hfill & \hfill \mathbf{O}\hfill & & & & & & & & & & & & & & & & & & \\ \hfill \mathbf{O}\hfill & \hfill \mathbf{O}\hfill & & & & & & & & & & & & & & & & & & \end{array}\right]{\mathit{\text{Q}}}^{\mathit{\text{t}}\mathit{\text{t}}},$$

(46)

with subsets ${\mathit{\text{s}}}^{\prime }\subset \mathit{\text{s}}$ and ${\mathit{\text{t}}}^{\prime }\subset \mathit{\text{t}}$. Here, ${\mathit{\text{P}}}^{\mathit{\text{s}}\mathit{\text{s}}}$ and ${\mathit{\text{Q}}}^{\mathit{\text{t}}\mathit{\text{t}}}$ are permutation matrices, ${\mathit{\text{Z}}}_{\text{E}}^{{\mathit{\text{s}}}^{\prime }{\mathit{\text{t}}}^{\prime }}$ denotes the matrix subblock corresponding to nonzero interactions. Note that, depending on the relative position between the two cubes considered, there can be many zero interactions in ${\mathit{\text{Z}}}_{\text{E}}^{\mathit{\text{s}}\mathit{\text{t}}}$, due to the sparsity pattern of ${\mathit{\text{Z}}}_{\text{E}}$. Therefore, the matrix compression is in fact applied to ${\mathit{\text{Z}}}_{\text{E}}^{{\mathit{\text{s}}}^{\prime }{\mathit{\text{t}}}^{\prime }}$ rather than ${\mathit{\text{Z}}}_{\text{E}}^{\mathit{\text{s}}\mathit{\text{t}}}$, and we have

$${\mathit{\text{Z}}}_{\text{E}}^{{\mathit{\text{s}}}^{\prime }{\mathit{\text{t}}}^{\prime }}={\mathit{\text{U}}}_{\text{E}}{\mathit{\text{V}}}_{\text{E}},$$

(47)

where $\left|{\mathit{\text{s}}}^{\prime }\right|$ and $\left|{\mathit{\text{t}}}^{\prime }\right|$ denote the sizes of index sets ${\mathit{\text{s}}}^{\prime }$ and ${\mathit{\text{t}}}^{\prime }$, respectively. Meanwhile, the dimension of the matrix ${\mathit{\text{U}}}_{\text{E}}$ is $\left|{\mathit{\text{s}}}^{\prime }\right|\times \mathrm{\Gamma }$, and the dimension of the matrix ${\mathit{\text{V}}}_{\text{E}}$ is $\mathrm{\Gamma }\times \left|{\mathit{\text{t}}}^{\prime }\right|$. Here, $\mathrm{\Gamma }$ denotes the compression rank of ${\mathit{\text{Z}}}_{\text{E}}^{{\mathit{\text{s}}}^{\prime }{\mathit{\text{t}}}^{\prime }}$ for a given tolerance. Ultimately, the compressed matrices ${\mathit{\text{U}}}_{\text{E}}$ and ${\mathit{\text{V}}}_{\text{E}}$ in (47) are stored, and are applied to the surface current vector I in an efficient manner. For many practical problems that involve multiscale discretization, the matrix block ${\mathit{\text{Z}}}_{\text{E}}^{{\mathit{\text{s}}}^{\prime }{\mathit{\text{t}}}^{\prime }}$ can be effectively compressed, and the memory consumption for ${\mathit{\text{U}}}_{\text{E}}$ and ${\mathit{\text{V}}}_{\text{E}}$ are greatly reduced compared to that for ${\mathit{\text{Z}}}_{\text{E}}^{\mathit{\text{s}}\mathit{\text{t}}}$.

III. NUMERICAL RESULTS

In this section, the computational performance of the proposed hybrid algorithm is demonstrated through several numerical examples. The conjugate gradient method (CG) [35] is used as the iterative solver, and the tolerance is set to 1.0E-3. The diagonal preconditioner is used to improve the convergence. Unless otherwise stated, the grid spacings (i.e., $h_x$, $h_y$, $h_z$) for the FFT acceleration are by default set to the common value $h=$ 0.1 $\lambda$, and the expansion orders (i.e., $M_x$, $M_y$, $M_z$) are by default set to the common value $M=3$. Furthermore, the well-known FFTW [34] is used to perform the FFT. For clarity, in the following, the traditional pre-splitting-FFT-based algorithm which employs the uncompressed ${\mathit{\text{Z}}}_{\text{E}}$ will be termed as the PSG-FFT algorithm, and the proposed hybrid fast algorithm (enhanced with the ACA) will be referred to as the PSG-FFT-ACA algorithm.

A. Sphere

Figure 7: Mesh configurations of the PEC sphere. The radius of the sphere is 1.0 $\lambda$. Here, about a quarter of the spherical surface is discretized with a relatively dense mesh of 0.02 $\lambda$, and the rest of the surface is discretized with a regular mesh of 0.1 $\lambda$. For clarity, the enlarged view of the multiscale mesh is also visualized in the above.

A PEC sphere of radius 1.0 $\lambda$ is first considered. The surface of the sphere is discretized with multiscale triangular meshes, resulting in 34,764 RWG functions. In particular, as shown in Fig. 7, about a quarter of the surface is discretized with dense meshes of edge lengths about 0.02 $\lambda$, and the rest is discretized with regular meshes of edge lengths about 0.1 $\lambda$.

For this example, the standard MoM (without using any acceleration technique), the traditional PSG-FFT algorithm, and the proposed PSG-FFT-ACA algorithm are adopted as the the solution methods. After solving the corresponding matrix equations as formulated in (4), the surface current coefficients (i.e, the solutions of equation (4)) ${\mathit{\text{I}}}^{\text{MoM}}$, ${\mathit{\text{I}}}^{\text{PSG-FFT}}$, and ${\mathit{\text{I}}}^{\text{PSG-FFT-ACA}}$ are then obtained, respectively. Here, it should be noticed that both ${\mathbf{I}}^{\text{PSG-FFT}}$ and ${\mathit{\text{I}}}^{\text{PSG-FFT-ACA}}$ depend on the choice of the pre-splitting threshold $\delta$. However, when the standard MoM is used, $\delta$ is not involved, and ${\mathit{\text{I}}}^{\text{MoM}}$ is therefore independent of $\delta$. Thus, the solution ${\mathit{\text{I}}}^{\text{MoM}}$ can be used as a proper reference, to evaluate the accuracy of ${\mathit{\text{I}}}^{\text{PSG-FFT}}$ and ${\mathit{\text{I}}}^{\text{PSG-FFT-ACA}}$ for varying $\delta$. The relative errors of the surface current coefficients are defined as follows

$$ℰ=\frac{{\parallel {\mathit{\text{I}}}^{\text{Test}}-{\mathit{\text{I}}}^{\text{Reference}}\parallel }_2}{{\parallel {\mathit{\text{I}}}^{\text{Reference}}\parallel }_2},$$

(48)

where ${\mathit{\text{I}}}^{\text{Test}}$ denotes the surface current coefficients to be tested, ${\mathit{\text{I}}}^{\text{Reference}}$ denotes some specified surface current coefficients to be used as the reference, and $\parallel \cdot {\parallel }_2$ denotes the 2-norm of vector. Based on (48), the relative errors of both the traditional PSG-FFT algorithm and the proposed PSG-FFT-ACA algorithm for several different $\delta$ are calculated and shown in Fig. 8. We can see that both algorithms can achieve accurate results and show similar accuracy level. In particular, within the range of the $\delta$ between 0.15 $\lambda$ and 0.65 $\lambda$, the errors of both algorithms decrease with increasing $\delta$. Furthermore, the accuracies of both algorithms are not strongly sensitive with respect to the $\delta$ considered.

Figure 8: Relative errors $ℰ$ of the surface currents coefficients I with respect to different splitting threshold $\delta$. The relative errors of both the traditional PSG-FFT algorithm and the proposed PSG-FFT-ACA algorithm are illustrated. The surface currents coefficients obtained using the standard MoM are used as the reference.

In the following, the choice of $\delta$ is discussed. In principle, $\delta$ can be of any value. However, in practice, a moderate value of $\delta$ is suggested, due to the following considerations.

• If $\delta$ is chosen be relatively large, the number of nonzero elements of the evanescent matrix ${\mathit{\text{Z}}}_{\text{E}}$ will then become relatively large. Clearly, in this case, the corresponding computational cost will increase. In Table 1, the computational costs with respect to $\delta$ are recorded. It can be observed that the time and memory consumptions of both the traditional PSG-FFT algorithm and the proposed PSG-FFT-ACA algorithm increase as $\delta$ increase. Hence, due to efficiency considerations, $\delta$ is not suggested to be too large.

• If $\delta$ is chosen to be relatively small, the number of nonzero elements of the evanescent matrix ${\mathit{\text{Z}}}_{\text{E}}$ will become relatively small, and the computational cost is correspondingly reduced. However, from Fig. 8, it is observed that the accuracy of the solution is slightly reduced when $\delta$ becomes smaller. Hence, due to accuracy considerations, $\delta$ is not suggested to be too small as well.

Based on the above considerations, the $\delta$ between 0.15 $\lambda$ and 0.65 $\lambda$ can be proper choice. In practice, by tuning the splitting threshold $\delta$ within this range, we can trade off between accuracy and computational cost.

Table 1: Time and memory for the sphere model with different splitting threshold $\delta$

$\delta (\mathit{\lambda })$	PSG-FFT		PSG-FFT-ACA
		Memory (GB)	CPU Time $(\mathbf{m})$	Memory (GB)	CPU Time $(\mathbf{m})$
$0.15(\lambda )$	0.9924	0.6601	0.9917	0.7087
$0.25(\lambda )$	1.9707	0.8078	1.6581	0.8167
$0.35(\lambda )$	3.1454	1.1445	2.2623	0.8679
$0.45(\lambda )$	4.6662	1.2814	2.9167	0.9177
$0.55(\lambda )$	6.3642	2.0607	3.5851	0.9812
$0.65(\lambda )$	8.3831	1.8004	4.2083	1.0390
$0.75(\lambda )$	10.514	2.1802	4.3177	1.0939
$0.85(\lambda )$	11.528	2.4548	5.1668	1.2197

Besides, the well-studied IE-FFT algorithm [17], being a typical instance of the pre-correction-based methods, is also used to calculate this example. Similar to the configurations of the pre-splitting-based methods employed, the grid spacing $h$ for the IE-FFT is set to 0.1 $\lambda$ and the expansion order $M$ is set to 3. Then, using ${\mathit{\text{I}}}^{\text{MoM}}$ as the reference, the relative error of the obtained surface current coefficients ${\mathit{\text{I}}}^{\text{IE-FFT}}$ is calculated, which is 2.96E-3. Referring to Fig. 8, we can see that such error is close to the error of the PSG-FFT(-ACA) algorithm with $\delta$ about 0.45 $\lambda$. Moreover, for the IE-FFT algorithm, the memory requirement is $3.09718$ G and the CPU time required is $1.08727$ m. Referring to Table 1, we can see that such cost is slightly lower than that of the PSG-FFT algorithm with $\delta$ about 0.45 $\lambda$. From the above numerical results, it can be known that the computational performance of the IE-FFT algorithm during the solution stage is slightly better than that of the PSG-FFT algorithm. However, it should be emphasized that such slight advantage of the IE-FFT algorithm is at the price of the additional grid-based numerical correction operations during the setup (precomputation) stage.

B. Monopole antenna on top of large platform

In the following, a more realistic example is considered, i.e., a monopole antenna on top of a large platform, as shown in Fig. 9. Clearly, this is a typical multiscale structure and often encountered in many practical scenarios. The width and height of the platform are 4.0 $\lambda$ and 0.3 $\lambda$, respectively. Moreover, the width and height of the monopole antenna are 0.06 $\lambda$ and 0.6 $\lambda$, respectively. The object surface is discretized with the multiscale triangular meshes. Specifically, the monopole antenna is discretized with meshes of edge length about 0.01 $\lambda$, while the large platform is discretized with meshes of edge length about 0.1 $\lambda$. With such mesh, a total of 24,807 RWG functions are generated.

Figure 9: Mesh configurations of a small monopole antenna on top of a large platform. A multiscale mesh is used to discretize this object. The monopole antenna is discretized using dense mesh of about 0.01 $\lambda$, while the large platform is discretized with regular mesh ofabout 0.1 $\lambda$.

Figure 10: Relative errors $ℰ$ of the surface currents coefficients with respect to different grid spacings $h$ and expansion orders $M$. The relative errors of the proposed PSG-FFT-ACA algorithm are illustrated. The surface currents coefficients obtained using the standard MoM are used as the reference.

The standard MoM and the proposed PSG-FFT-ACA algorithm are employed to calculate this example. Here, the pre-splitting threshold $\delta$ is chosen to be 0.35 $\lambda$. Using ${\mathit{\text{I}}}^{\text{MoM}}$ as the reference, the relative errors of ${\mathit{\text{I}}}^{\text{PSG-FFT-ACA}}$ are calculated for different grid spacings $h$ and different expansion orders $M$. For clarity, the relative errors are illustrated in Fig. 10. The corresponding time and memory consumptions are recorded in Table 2. It can be clearly seen that, for decreasing grid spacing $h$ and increasing expansion order $M$, the relative error decreases and the corresponding computational cost increases. In other words, by tuning $h$ and $M$, we can trade off between accuracy and computational cost.

Table 2: Time and memory for the monopole-platform model with different grid spacing $h$ and expansion order $M$

$h(\lambda )$	$M=2$		$M=3$
		Memory (GB)	CPU Time $(\mathbf{m})$	Memory (GB)	CPU Time $(\mathbf{m})$
$0.05(\lambda )$	1.0074	0.7863	1.1094	0.8256
$0.10(\lambda )$	0.8643	0.2280	0.9662	0.2792
$0.15(\lambda )$	0.8311	0.1554	0.9304	0.2106
$0.20(\lambda )$	0.8243	0.1155	0.9261	0.1801

C. Cone-shape object

To further study the effectiveness of the matrix compression enhancement, a PEC cone model is then considered, as shown in Fig. 11. The height of the cone is 3.0 $\lambda$. The radii of the top and the bottom faces are 0.03 $\lambda$ and 0.6 $\lambda$, respectively. The surface is discretized with a gradually-varied meshes. In particular, three mesh configurations are considered. For the three cases, the edge length for the regular part of the mesh is fixed to 0.1 $\lambda$, while the edge lengths for the dense part of the mesh are set to 0.01 $\lambda$, 0.005 $\lambda$, and 0.0025 $\lambda$, respectively. For clarity, the mesh settings considered are illustrated in Fig. 11.

Figure 11: Mesh configurations of the PEC cone. A gradually-varied mesh is used to discretize the cone-shape object. The edge length for the regular part of the mesh is 0.1 $\lambda$. The edge lengths for the dense part of the mesh are 0.01 $\lambda$, 0.005 $\lambda$, and 0.0025 $\lambda$, respectively. For clarity, the enlarged views are also shown in the above.

Table 3: Time and memory for the cone model of case 1

Algorithm	Memory (GB)		CPU Time
		$Z_{\mathrm{E}}$	Total	$Z_{\mathrm{E}}\cdot I$ $(\mathrm{s})$	Total $(\mathbf{m})$
PSG-FFT	0.3313	0.6286	12.128	0.6955
PSG-FFT-ACA	0.2418	0.4738	8.3050	0.6194

Table 4: Time and memory for the cone model of case 2

Algorithm	Memory (GB)		CPU Time
		$Z_{\mathrm{E}}$	Total	$Z_{\mathrm{E}}\cdot I$ $(\mathrm{s})$	Total $(\mathrm{m})$
PSG-FFT	0.9611	1.6176	69.751	2.3978
PSG-FFT-ACA	0.4935	0.8778	35.127	1.8475

For each of these mesh configurations, we solve the corresponding equations first using the PSG-FFT algorithm which employs the uncompressed ${\mathit{\text{Z}}}_{\text{E}}$. Then, we recalculate these problems using the PSG-FFT-ACA algorithm which compresses ${\mathit{\text{Z}}}_{\text{E}}$ with the ACA. Here, the splitting threshold $\delta$ is chosen to be 0.35 $\lambda$. Notice that, with such configuration of $\delta$, all the interactions corresponding to ${\mathit{\text{Z}}}_{\text{E}}$ occur within the subwavelength regime. Thus, the ACA compression can be applied effectively [21, 37]. In Table  3, 4, 5, we tabulate the time and memory consumptions of both the PSG-FFT algorithm and the PSG-FFT-ACA algorithm for the three mesh cases considered. It can be seen that, for all the cases, the PSG-FFT-ACA algorithm shows improved computational performances relative to the PSG-FFT algorithm. Further inspections show that the degree of improvements becomes more prominent when the dense mesh is finer. To make this more clear, we calculate the compression ratios of the evanescent matrix ${\mathit{\text{Z}}}_{\text{E}}$ and use this as an indicator to measure the degree of improvements of the matrix compression enhancement. In Fig. 12, the compression ratios with respect to different mesh configurations are illustrated. We can clearly see that the compression ratio becomes higher when the dense mesh becomes finer. That is, the effectiveness of the proposed hybrid algorithm (with matrix compression enhancement) becomes more prominent when the ratio between the regular mesh size and the dense mesh size is larger.

Table 5: Time and memory for the cone model of case 3

Algorithm	Memory (GB)		CPU Time
		$Z_{\mathrm{E}}$	Total	$Z_{\mathrm{E}}\cdot \mathit{I}$ $(\mathbf{m})$	Total $(\mathbf{m})$
PSG-FFT	2.6465	4.2615	181.47	4.4874
PSG-FFT-ACA	1.1579	2.0570	78.678	2.8009

Figure 12: Compression ratios of the cone with different mesh configurations.

D. Aircraft model

A relatively large PEC aircraft model with multiscale meshes is finally considered, as shown in Fig. 13. In this model, both the length and the wingspan of the aircraft are about 16 $\lambda$. The surface is discretized with triangular meshes, resulting in 120,072 RWG functions. Here, a part of the aircraft surface on the wings is discretized with dense meshes of edge lengths about 0.02 $\lambda$, and the rest of the aircraft surface is discretized with regular meshes of edge lengths about 0.1 $\lambda$.

The bistatic RCS curves obtained by the PSG-FFT algorithm and the PSG-FFT-ACA algorithm are illustrated in Fig. 14. The pre-splitting threshold $\delta$ is here chosen to be 0.35 $\lambda$. The bistatic RCS curve obtained by the well-tested IE-FFT algorithm [17] is also shown in the figure and used as the reference to calculate the root mean square error (RMSE). We can see that these curves agree well with each other.

Figure 13: Geometry model and corresponding multiscale mesh of the PEC aircraft. Here, a part of the aircraft surface on the wings are discretized with a relatively dense mesh of 0.02 $\lambda$, while the rest of the aircraft surface is discretized with a regular mesh of 0.1 $\lambda$. The enlarged view of the multiscale mesh is also highlighted in the above.

Figure 14: Bistatic RCS curves of the PEC aircraft model. The result obtained by the IE-FFT is used as the reference. The RMSE of the PSG-FFT and the PSG-FFT-ACA are calculated.

The time and memory consumptions of both the PSG-FFT algorithm and the PSG-FFT-ACA algorithm are tabulated in Table  6. From the table, we can see that the memory requirement of the evanescent matrix ${\mathit{\text{Z}}}_{\text{E}}$ in the PSG-FFT-ACA algorithm is reduced by about $34.1218\%$, compared with that of the PSG-FFT algorithm. Meanwhile, the time cost for performing ${\mathit{\text{Z}}}_{\text{E}}\cdot \mathit{\text{I}}$ in the whole solution process is reduced by about $37.52\%$. Thus, the time and memory cost related to the evanescent matrix ${\mathit{\text{Z}}}_{\text{E}}$ can be greatly reduced, as a result of the matrix compression enhancement.

Nevertheless, the total cost of the overall algorithm are determined not only by the cost related to the evanescent matrix ${\mathit{\text{Z}}}_{\text{E}}$, but also determined by the cost related to the propagating matrix ${\mathit{\text{Z}}}_{\text{P}}$. For this relatively large object, the dense mesh region only occupies a relatively small portion of the overall object surface, and the ratio between the edge length (0.1 $\lambda$) of the regular mesh and the edge length (0.02 $\lambda$) of the dense mesh is here not extremely large. Thus, the time cost for performing ${\mathit{\text{Z}}}_{\text{P}}\cdot \mathit{\text{I}}$ still dominates the time cost for performing $\mathit{\text{Z}}\cdot \mathit{\text{I}}$, and the reduction of the time cost for performing ${\mathit{\text{Z}}}_{\text{E}}\cdot \mathit{\text{I}}$ therefore only has a minor influence on the total time cost for performing $\mathit{\text{Z}}\cdot \mathit{\text{I}}$. As a result, although the total time can be reduced due to the reduced cost related to ${\mathit{\text{Z}}}_{\text{E}}$, the reduction of the total time does not look very striking for this example.

Table 6: Time and memory for the aircraft model

Algorithm	Memory (GB)		CPU Time
		$Z_{\mathrm{E}}$	Total	$Z_{\mathrm{E}}\cdot I$ $(\mathrm{s})$	Total $(\mathbf{m})$
PSG-FFT	3.5558	7.7408	76.910	23.580
PSG-FFT-ACA	2.3425	5.4033	48.053	21.431

IV. CONCLUSIONS

In this paper, a hybrid fast algorithm has been established for efficiently analyzing multiscale problems. In this algorithm, the Green’s function is a priori split into a singular but spatially localized part, and a smooth, oscillatory but bandlimited part. Then, the fundamental blocks for acceleration, i.e., the FFT and the matrix compression method (ACA), are applied to these two parts, respectively. Compared with the traditional methods which only employ the FFT for acceleration, the proposed hybrid algorithm can maintain low memory consumption in multiscale cases without compromising time cost, thus manifesting itself as a favorable multiscale extension of the traditional pure FFT-based methods. Furthermore, the required auxiliary geometric data structures herein can be constructed independent of each other, thus permitting the two kinds of acceleration methods to be connected seamlessly in a concise and elegant manner.

On the basis of the present work, several future directions may be explored. First, except for the Lagrange interpolation used in this work, other more accurate interpolation schemes (such as the Gaussian interpolation [15] or the fitting techniques [12]) can be employed to further improve the computational performance of the overall algorithm. Moreover, although only the PEC problems are considered here, extensions to the dielectric problems (especially inhomogeneous cases [16]) are clearly feasible. Furthermore, effective parallelization (especially on heterogeneous architectures [13]) of the proposed hybrid algorithm can be interesting subject of future research.

APPENDIX

The derivative functions of $\mathrm{Re}\left[\mathrm{G}\left(R\right)\right]$, which are used in (18) and (19), are given below

$$\frac{d}{dR}\mathrm{Re}\left[\mathrm{G}\left(R\right)\right]=-\frac{\cos \left(kR\right)+kR\sin \left(kR\right)}{4\pi R^2},$$

(49)

$$\begin{array}{c}\hfill \frac{d^2}{d{R}^2}\mathrm{Re}\left[\mathrm{G}\left(R\right)\right]=\hfill \\ \hfill -\frac{\left(k^2{R}^2-2\right)\cos \left(kR\right)-2kR\sin \left(kR\right)}{4\pi R^3}.\hfill \end{array}$$

(50)

The derivative functions of $\mathrm{Re}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(R\right)\right]$, which are used in (30) and (31), are given below

$$\begin{array}{c}\hfill \frac{d}{dR}\mathrm{Re}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(R\right)\right]=\hfill \\ \hfill -\frac{\left(k^2{R}^2-2\right)\cos \left(kR\right)-2kR\sin \left(kR\right)}{4\pi R^3},\hfill \end{array}$$

(51)

$$\begin{array}{c}\hfill \frac{d^2}{d{R}^2}\mathrm{Re}\left[\stackrel{\dot{}\dot{}\dot{}}{\mathrm{G}}\left(R\right)\right]=\hfill \\ \hfill -\frac{3\left(k^2{R}^2-2\right)\cos \left(kR\right)-kR\left(k^2{R}^2-6\right)\sin \left(kR\right)}{4\pi R^4}.\hfill \end{array}$$

(52)

ACKNOWLEDGMENT

This work was supported in part by the Fundamental Research Funds for the Central Universities (No. 2242022k30003), in part by the National Key Research and Development Project (No. 2019YFB1803202), and in part by the National Science Foundation of China (Nos. 62293492, 62131008, 61975177, U20A20164).

REFERENCES

[1] W. C. Chew and L. J. Jiang, “Overview of large-scale computing: the past, the present, and the future,” Proceedings of the IEEE, vol. 101, no. 2, pp. 227-241, Feb. 2013.

[2] M.-K. Li and W. C. Chew, “Multiscale simulation of complex structures using equivalence principle algorithm with high-order field point sampling scheme,” IEEE Trans. Antennas Propag., vol. 56, no. 8, pp. 2389-2397, 2007.

[3] Z.-G. Qian and W. C. Chew, “Fast full-wave surface integral equation solver for multiscale structure modeling,” IEEE Trans. Antennas Propag., vol. 57, no. 11, pp. 3594-3601, Nov. 2009.

[4] J. Zhu, S. Omar, and D. Jiao, “Solution of the electric field integral equation when it breaks down,” IEEE Trans. Antennas Propag., vol. 62, no. 8, pp. 4122-4134, Aug. 2014.

[5] R. F. Harrington, Field Computation by Moment Methods. New York: The MacMillian Co., 1968.

[6] J. M. Song and W. C. Chew, “Multilevel fast-multipole algorithm for solving combined field integral equations of electromagnetic scattering,” Microw. Opt. Technol. Lett., vol. 10, no. 1, pp. 14-19, Sep. 1995.

[7] J.-M. Song, C.-C. Lu, and W. C. Chew, “Multilevel fast multipole algorithm for electromagnetic scattering by large complex objects,” IEEE Trans. Antennas Propag., vol. 45, no. 10, pp. 1488-1493, Oct. 1997.

[8] W. C. Chew, J.-M. Jin, E. Michielssen, and J. Song, Fast and Efficient Algorithms in Computational Electromagnetics. Boston, MA: Artech House, 2001.

[9] M. Bleszynski, E. Bleszynski, and T. Jaroszewicz, “AIM: Adaptive integral method for solving large-scale electromagnetic scattering and radiation problems,” Radio Sci., vol. 31, no. 5, pp. 1225-1251, 1996.

[10] J. R. Phillips and J. K. White, “A precorrected-FFT method for electrostatic analysis of complicated 3-D structures,” IEEE Trans. Computer-Aided Design, vol. 16, no. 10, pp. 1059-1072, Oct. 1997.

[11] S. M. Seo and J.-F. Lee, “A fast IE-FFT algorithm for solving PEC scattering problems,” IEEE Trans. Magn., vol. 41, no. 5, pp. 1476-1479, May2005.

[12] J.-Y. Xie, H.-X. Zhou, W. Hong, W.-D. Li, and G. Hua, “A highly accurate FGG-FG-FFT for the combined field integral equation,” IEEE Trans. Antennas Propag., vol. 61, no. 9, pp. 4641-4652, Sep. 2013.

[13] S. Peng and C.-F. Wang, “Precorrected-FFT method on graphics processing units,” IEEE Trans. Antennas Propag., vol. 61, no. 4, pp. 2099-2107, Apr. 2013.

[14] L.-W. Li, Y.-J. Wang, and E.-P. Li, “MPI-based parallelized precorrected FFT algorithm for analyzing scattering by arbitrarily shaped three-dimensional objects,” Electromagn. Waves Appl., vol. 17, pp. 1489-1491, Jan. 2003.

[15] B. Lai, X. An, H.-B. Yuan, Z.-S. Chen, C.-M. Huang, and C.-H. Liang, “A novel Gaussian interpolation formula-based IE-FFT algorithm for solving EMscattering problems,” Microw. Opt. Tech. Lett., vol. 51, no. 9, pp. 2233-2236, Sep.2009.

[16] X.-C. Nie, L.-W. Li, N. Yuan, T. S. Yeo, and Y.-B. Gan, “Precorrected-FFT solution of the volume integral equation for 3-D inhomogeneous dielectric objects,” IEEE Trans. Antennas Propag., vol. 53, no. 1, pp. 313-320, Jan. 2005.

[17] J.-Y. Xie, H.-X. Zhou, W.-D. Li, and W. Hong, “IE-FFT for the combined field integral equation applied to electrically large objects,” Microw. Opt. Tech. Lett., vol. 54, no. 4, pp. 891-896, Apr.2012.

[18] S. Peng, C.-F. Wang, and J.-F. Lee, “Analyzing PEC scattering structure using an IE-FFT algorithm,” ACES Journal, vol. 24, no. 2, Apr. 2009.

[19] L. J. Jiang and W. C. Chew, “A mixed-form fast multipole algorithm,” IEEE Trans. Antennas Propag., vol. 53, no. 12, pp. 4145-4156, Dec.2005.

[20] X.-M. Pan, J.-G. Wei, Z. Peng, and X.-Q. Sheng, “A fast algorithm for multiscale electromagnetic problems using interpolative decomposition and multilevel fast multipole algorithm,” Radio Sci., vol. 47, no. 1, Feb. 2012.

[21] L.-F. Ma, Z. Nie, J. Hu, and S.-Q. He, “Combined MLFMA-ACA algorithm application to scattering problems with complex and fine structure,” Asia Pacific Microwave Conference, pp. 802-805,2009.

[22] W.-B. Kong, H.-X. Zhou, K.-L. Zheng, and W. Hong, “Analysis of multiscale problems using the MLFMA with the assistance of the FFT-based method,” IEEE Transactions on Antennas and Propagation, vol. 63, no. 9, pp. 4184-4188,2015.

[23] J. S. Zhao and W. C. Chew, “Three dimensional multilevel fast multipole algorithm from static to electrodynamic,” Microw. Opt. Technol. Lett., vol. 26, no. 1, pp. 43-48, 2000.

[24] L. J. Jiang and W. C. Chew, “Low-frequency fast inhomogeneous plane-wave algorithm,” Microw. Opt. Tech. Lett., vol. 40, no. 2, pp. 117-122, Jan. 2004.

[25] T. Xia, L. L. Meng, Q. S. Liu, H. H. Gan, and W. C. Chew, “A low-frequency stable broadband multilevel fast multipole algorithm using plane wave multipole hybridization,” IEEE Trans. Antennas Propag., vol. 66, no. 11, pp. 6137-6145, Nov.2018.

[26] H. Wallen and J. Sarvas, “Translation procedures for broadband MLFMA,” Prog. Electromagn. Res., vol. 55, pp. 47-78, 2005.

[27] D. Wulf and R. Bunger, “An efficient implementation of the combined wideband MLFMA/LF-FIPWA,” IEEE Trans. Antennas Propag., vol. 57, no. 2, pp. 467-474, Feb. 2009.

[28] D. T. Schobert and T. F. Eibert, “Fast integral equation solution by multilevel Green’s function interpolation combined with multilevel fast multipole method,” IEEE Trans. Antennas Propag., vol. 60, no. 9, pp. 4458-4463, Sep. 2012.

[29] I. Bogaert, J. Peeters, and F. Olyslager, “A nondirective plane wave MLFMA stable at low frequencies,” IEEE Trans. Antennas Propag., vol. 56, no. 12, pp. 3752-3767, Dec. 2008.

[30] O. Ergul and B. Karaosmanoglu, “Broadband multilevel fast multipole algorithm based on an approximate diagonalization of the Green’s function,” IEEE Trans. Antennas Propag., vol. 63, no. 7, pp. 3035-3041, July 2015.

[31] W.-B. Kong, H.-X. Zhou, K.-L. Zheng, X. Mu, and W. Hong, “FFT-based method with near-matrix compression,” IEEE Transactions on Antennas and Propagation, vol. 65, no. 11, pp. 5975-5983,2017.

[32] A. Bespalov, “Cost-effective solution of the boundary integral equations for 3D Maxwell problems,” Russ. J. Numer. Anal. Math. Modelling, vol. 14, no. 5, pp. 403-428, 1999.

[33] A. Bespalov, “On the use of a regular grid for implementation of boundary integral methods for wave problems,” Russ. J. Numer. Anal. Math. Modelling, vol. 15, no. 6, pp. 469-488. 2000.

[34] M. Frigo and S. Johnson, FFTW Manual [Online]. Available: http://www.fftw.org/

[35] R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. M. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, and H. Van der Vorst, Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. Philadelphia, PA, USA: SIAM,1994.

[36] M. Bebendorf and S. Rjasanow, “Adaptive low-rank approximation of collocation matrices,” Computing, vol. 70, no. 1, pp. 1-24, 2003.

[37] K. Zhao, M. N. Vouvakis, and J.-F. Lee, “The adaptive cross approximation algorithm for accelerated method of moments computations of EMC problems,” IEEE Trans. Electromag. Compat., vol. 47, no. 4, pp. 763-773, Nov. 2005.

[38] H.-X. Zhou, G.-Y. Zhu, W.-B. Kong, and W. Hong, “An upgraded ACA algorithm in complex field and its statistical analysis,” IEEE Trans. Antennas Propag., vol. 65, no. 5, pp. 2734-2739, May2017.

[39] S. M. Rao, D. R. Wilton, and A. W. Glisson, “Elettromagnetic scattering by surfaces of arbitrary shape,” IEEE Trans. Antennas Propag., vol. 30, no. 3, pp. 409-418, May 1982.

BIOGRAPHIES

Guang-Yu Zhu received the B.D. in electrical engineering from Nanjing Agricultural University, Nanjing, China, in 2012, and the Ph.D. degree in radio engineering from Southeast University, Nanjing, China, in 2020, respectively. He is currently a postdoctoral research fellow with the Department of Information Science and Electronic Engineering, Zhejiang University, Hangzhou, China. His research interests are in computational electromagnetics with focus on fast algorithms and domain decomposition methods.

Wei-Dong Li received the M.S. degree in mathematics and the Ph.D. degree in radio engineering from Southeast University, Nanjing, China, in 2003 and 2007, respectively. From January 2008 to January 2009, he was a visiting scholar with the Technische University Darmstadt, Germany. He is currently an associate professor with the State Key Laboratory of Millimeter Waves and the School of Information Science and Engineering, Southeast University. He has authored or coauthored over 40 technical articles. He serves as a reviewer for the IEEE Transactions on Antennas and Propagation and IET Microwave, Antennas and Propagation. His research interests are in computational EM with focus on integral-equation-based overlapped domain decomposition method, multilevel fast multipole algorithm, fast and accurate inter/extrapolation techniques, and basic functions.

Wei E.I. Sha received the B.S. and Ph.D. degrees in electronic engineering from Anhui University, Hefei, China, in 2003 and 2008, respectively. From July 2008 to July 2017, he was a post-doctoral research fellow and then a research assistant professor with the Department of Electrical and Electronic Engineering, University of Hong Kong, Hong Kong. From March 2018 to March 2019, he worked as a Marie Skodowska-Curie Individual Fellow with University College London, London, U.K. From October 2017, he joined the College of Information Science and Electronic Engineering, Zhejiang University, Hangzhou, China, where he is currently a tenure-tracked assistant professor. He has authored or coauthored 180 refereed journal articles, 150 conference publications (including five keynote talks and one short course), nine book chapters, and two books. His Google Scholar citation is 8193 with H-index of 45. His research interests include theoretical and computational research in electromagnetics and optics, focusing on the multiphysics and interdisciplinary research. His research involves fundamental and applied aspects in computational and applied electromagnetics, nonlinear and quantum electromagnetics, micro- and nano-optics, optoelectronic device simulation, and multiphysics modeling. Dr. Sha is a member of Optical Society of America (OSA). He served as a reviewer for 60 technical journals and technical program committee member for ten IEEE conferences. He was a recipient of the Applied Computational Electromagnetics Society (ACES) Technical Achievement Award in 2022 and PIERS Young Scientist Award in 2021. In 2015, he was awarded Second Prize of Science and Technology from Anhui Province Government, China. In 2007, he was awarded the Thousand Talents Program for Distinguished Young Scholars of China. He also received six Best Student Paper Prizes and one Young Scientist Award with his students. He also served as an associate editor for IEEE Journal on Multiscale and Multiphysics Computational Techniques, IEEE Open Journal of Antennas and Propagation, and IEEE Access.

Hou-Xing Zhou received the M.S. degree in mathematics from Southwest Normal University, Chongqing, China, in 1995, and the Ph.D. degree in radio engineering from Southeast University, Nanjing, China, in 2002. Since 2002, he has been with the State Key Laboratory of Millimeter Waves, Southeast University. He is currently a professor with the School of Information Science and Engineering, Southeast University. He has authored or coauthored 40 journal articles. His main research interests include numerical algorithms in computational EM, including the fast algorithm for spatial domain dyadic Green’s functions of stratified media, the multilevel fast multiple algorithm, the FFT-based fast algorithm, the IE-based domain decomposition method, the FEM-BI-based domain decomposition method, the higher order method of moments, and the parallel computation based on GPU/multicore-CPU platform. In addition, the full-wave simulation technology for the electromagnetic field distribution around the base-station is one of his current interests.

Wei Hong received the B.S. degree from the University of Information Engineering, Zhengzhou, China, in 1982, and the M.S. and Ph.D. degrees from Southeast University, Nanjing, China, in 1985 and 1988, respectively, all in radio engineering. Since 1988, he has been with the State Key Laboratory of Millimeter Waves and serves for the director of the lab, since 2003. In 1993, 1995, 1996, 1997, and 1998, he was a short-term visiting scholar with the University of California at Berkeley and at Santa Cruz, respectively. He is currently a professor with the School of Information Science and Engineering, Southeast University. He has been involved in numerical methods for electromagnetic problems, millimeter wave theory and technology, antennas, RF technology for wireless communications, and so on. He has authored and coauthored over 300 technical publications and two books. Dr. Hong was an elected IEEE MTT-S AdCom Member, from 2014 to 2016. He is a fellow of CIE, the vice president of the CIE Microwave Society and Antenna Society, and the chair of the IEEE MTT-S/AP-S/EMC-S Joint Nanjing Chapter. He was twice awarded the National Natural Prizes, thrice awarded the first-class Science and Technology Progress Prizes issued by the Ministry of Education of China and Jiangsu Province Government. He also received the Foundations for China Distinguished Young Investigators and for Innovation Group issued by NSF of China. He served as the associate editor for the IEEE Transactions on Microwave Theory and Techniques from 2007 to 2010, and one of the guest editors for the 5G special issue of the IEEE Transactions on Antennas and Propagation, in 2017.

ACES JOURNAL, Vol. 38, No. 9, 652–666
doi: 10.13052/2023.ACES.J.380904
© 2023 River Publishers