Efficient Improved Local Time Stepping with the Leapfrog Scheme for Transient 3-D Electromagnetic Analyses

Minxuan Li${}^{\mathbf{\text{1,2}}}$, Qingkai Wu${}^{\mathbf{\text{1}}}$, Zhongchao Lin${}^{\mathbf{\text{1}}}$, Yu Zhang${}^{\mathbf{\text{1}}}$ and Xunwang Zhao${}^{\mathbf{\text{1}}}$

¹School of Electrical Engineering
Xidian University, Xi’an, 710071, China
zclin@xidian.edu.cn
²China Electronic Product Reliability and Environmental Testing Institute
Guangzhou, 511300, China
lwdldj@163.com

Submitted On: June 29, 2022; Accepted On: September 7, 2023

ABSTRACT

An alternate boundary local time stepping (ABLTS) method is proposed for the discontinuous Galerkin time domain method for transient electromagnetic simulations to reduce the computation complexity of the local time stepping (LTS) method. The proposed method exhibites lower storage and time complexity than the conventional LF-LTS method . The stability, accuracy, and effectiveness of the ABLTS method are verified by applying it to the simulation of a resonator cavity and multi-layer microstrip antenna. The numerical results revealed that the developed method is effective for the transient electromagnetic simulation of antennas.

Index Terms: Antenna transient analysis, discontinuous Galerkin time domain method, high-order time integration, local time stepping.

I. INTRODUCTION

The discontinuous Galerkin time domain (DGTD) method is widely used in the transient simulation of antennas and microwave devices because of its high accuracy. The DGTD method introduces numerical flux into the finite element time domain (FETD) method [1–8] and exhibits higher accuracy and modeling flexibility than conventional time-domain methods such as finite difference time-domain (FDTD) method. The numerical flux decomposes the common basis function between neighboring elements, so that the governing equation can be established in one element and avoid building and decomposing large sparse matrices.

Runge-Kutta and leapfrog (LF) methods are the widely used explicit time integration schemes based on the DGTD method [6, 8, 9]. The LF scheme exhibits limited iteration steps. In the Courant–Friedrichs–Lewy (CFL) stability condition, the maximum time step size of explicit difference schemes is determined by the smallest element. However, unstructured meshes typically produce some distorted and small-size elements, which leads to a minimal global time step size [9]. Local time stepping (LTS) methods are typically used to address the multi-scale problem [10–12]. The LTS method classifies elements into many parts according to their size in space; fine and coarse regions have different time step sizes [8–14]. However, if a fine region is adjacent to a coarse region, the fine region requires the numerical flux of the adjacent coarse region, and the coarse region cannot provide fields at the exact time because of the large time step. The LTS method based on the arbitrary high order (ADER) scheme compensates for the errors of different time steps through the high-order Taylor expansion of time partial conductance in the coarse region. The ADER-LTS method is highly accurate but requires time synchronization when it provides field information as the output. Furthermore, frequent synchronization leads to complications in time integration steps, the ADER-LTS requires extra memory because of storing fields at multiple time steps. Another LTS method based on the leapfrog (LF) scheme is frequently discussed [8–12]. The LF scheme presents second-order time accuracy and simpler iterations than the ADER scheme. When adjacent element fields are not obtained at the exact time, the accuracy of the LF scheme degrades to first-order and some interpolation methods [8, 11], which increase iteration complexity and memory, are used to compensate for the accuracy. However the interpolation slows the solving speed and needs to meet a more stringent stability condition.

In this study, we propose an alternate boundary LTS (ABLTS) method based on the LF-LTS method reported in [8]. The ABLTS method is applied to the DGTD method with hierarchical vector basis functions, avoids interpolation between coarse and fine regions, and maintains the concise iterations of the LF scheme. The storage and time complexity of this method are lower than those of the interpolation LF-LTS method, which improves the universality of the LTS method. The simulation of a resonator cavity and patch antenna showed that the ABLTS method presents the high accuracy, fast speed, and low memory as the interpolation LF-LTS method.

II. PRINCIPLE AND FORMULATIONS

A. Numerical discretization of LTS-DGTD

The isotropic Maxwell’s curl equations in 3-D space without sources and lossless are

$$\{\begin{array}{cccccccccccccccccccc}\hfill \nabla \times \mathit{\text{H}}=\epsilon \frac{\partial \mathit{\text{E}}}{\partial t}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \nabla \times \mathit{\text{E}}=-\mu \frac{\partial \mathit{\text{H}}}{\partial t}\hfill & & & & & & & & & & & & & & & & & & & \end{array},$$

(1)

where $\epsilon$ is the permittivity, $\mu$ is the permeability. markedasmathbold italicE and markedasmathbold italicH are expended with the second-order hierarchical vector basis function reported in [6].

Using the Galerkin finite elements approach reported in [1], (1) can be expressed as the weak form with integration by parts as follows:

$$\{\begin{array}{cccccccccccccccccccc}\hfill {\int }_{{\mathrm{\Omega }}_i}\epsilon \frac{\partial \mathit{\text{E}}}{\partial t}\cdot {\mathit{\text{N}}}_i\mathit{d}\nu -{\int }_{{\mathrm{\Omega }}_i}\nabla \times {\mathit{\text{N}}}_i\cdot \mathit{\text{H}}\mathit{d}\nu =\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\int }_{{\mathrm{\Gamma }}_i}\left(\widehat{\mathit{\text{n}}}\times \mathit{\text{H}}\right)\cdot {\mathit{\text{N}}}_i\mathit{d}s,\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\int }_{{\mathrm{\Omega }}_i}\mu \frac{\partial \mathit{\text{H}}}{\partial t}\cdot {\mathit{\text{N}}}_i\mathit{d}\nu +{\int }_{{\mathrm{\Omega }}_i}\nabla \times {\mathit{\text{N}}}_i\cdot \mathit{\text{E}}\mathit{d}\nu =\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill -{\int }_{{\mathrm{\Gamma }}_i}\left(\widehat{\mathit{\text{n}}}\times \mathit{\text{E}}\right)\cdot {\mathit{\text{N}}}_i\mathit{d}s,\hfill & & & & & & & & & & & & & & & & & & & \end{array}$$

(2)

where $\mathrm{\Omega }$ is the computational domain divided into tetrahedrons ${\mathrm{\Omega }}_i$, ${\mathrm{\Gamma }}_i$ is the boundary of ${\mathrm{\Omega }}_i$, and $\widehat{\mathit{\text{n}}}$ is the normal vector located on ${\mathrm{\Gamma }}_i$ and points to the outside, ${\mathit{\text{N}}}_i$ is the weight functions of ${\mathrm{\Omega }}_i$. The DGTD method introduces numerical fluxes to evaluate the integration over tetrahedron interfaces. As a common numerical flux with excellent convergence, the upwind flux is introduced into (2), and the semi-discrete form with the upwind flux is expressed as follows [1]:

$$\begin{array}{cccccccccccccccccccc}& \hfill \epsilon {\mathbf{M}}_i\frac{\partial {\mathit{\text{E}}}_i}{\partial t}={\mathbf{S}}_i{\mathit{\text{H}}}_i+\hfill & & & & & & & & & & & & & & & & & & \\ & \hfill \sum _{p=1}^4(k_h\mathbf{F}{}_{{\mathrm{\Gamma }}_p}{}_i({\mathit{\text{H}}}_j-{\mathit{\text{H}}}_i)-v_e\mathbf{G}{}_{{\mathrm{\Gamma }}_p}{}_i({\mathit{\text{E}}}_j-{\mathit{\text{E}}}_i)),\hfill & & & & & & & & & & & & & & & & & & \\ & \hfill \mu {\mathbf{M}}_i\frac{\partial {\mathit{\text{H}}}_i}{\partial t}=-{\mathbf{S}}_i{\mathit{\text{E}}}_i-\hfill & & & & & & & & & & & & & & & & & & \\ & \hfill \sum _{p=1}^4(k_e\mathbf{F}{}_{{\mathrm{\Gamma }}_p}{}_i({\mathit{\text{E}}}_j-{\mathit{\text{E}}}_i)+v_h\mathbf{G}{}_{{\mathrm{\Gamma }}_p}{}_i({\mathit{\text{H}}}_j-{\mathit{\text{H}}}_i)),\hfill & & & & & & & & & & & & & & & & & & \end{array}$$

(3)

where $k_e$, $v_e$, $k_h$, and $v_h$ are the upwind flux coefficients [5], markedasmathboldM is the mass matrix, markedasmathboldS is the stiffness matrix, and ${\mathbf{F}}_{\mathrm{\Gamma }}$ and ${\mathbf{G}}_{\mathrm{\Gamma }}$ are the flux matrices.

(3) can be expressed in a highly concise form as follows:

$$\frac{\partial \mathit{\text{u}}}{\partial t}=\mathbf{A}\left(\mathit{\text{u}}\right),$$

(4)

where $\mathit{\text{u}}={[\mathit{\text{E}},\mathit{\text{H}}]}^T$, $\mathbf{A}$ denotes the terms of (3) other than the time partial term. According to the ADER scheme reported in [12], markedasmathbold italicu can be expanded into Taylor series as follows:

$$\mathit{\text{u}}\left(t+\mathrm{\Delta }t\right)=\sum _{p=0}^{N_t}\frac{\mathrm{\Delta }{t}^p}{p!}\frac{{\partial }^p\mathit{\text{u}}\left(t\right)}{\partial t^p},$$

(5)

where $N_t$ is the order in the ADER scheme. The first order ADER scheme is expressed as follows:

$$\frac{\partial \mathit{\text{u}}\left(t\right)}{\partial t}=\frac{\mathit{\text{u}}\left(t+\mathrm{\Delta }t\right)-\mathit{\text{u}}\left(t\right)}{\mathrm{\Delta }t}+O(t^2)=\mathbf{A}(\mathit{\text{u}}).$$

(6)

For isotropic media, $\mathbf{A}\left(\mathit{\text{u}}\right)=\mathbf{A}\mathit{\text{u}}(t)$, this scheme becomes an explicit format and can be derived into the LF scheme using (6) as follows:

$$\begin{array}{cccccccccccccccccccc}\hfill \frac{\mathit{\text{H}}_i{}^{n+\frac{1}{2}}-\mathit{\text{H}}_i{}^{n-\frac{1}{2}}}{\mathrm{\Delta }t}={\mu }^{-1}\mathbf{M}_i{}^{-1}\left(-{\mathbf{S}}_i\mathit{\text{E}}_i{}^n+f_h{}_i(\mathit{\text{u}})\right)\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill =\mathbf{A}_h{}_i\mathit{\text{u}}(t),\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \frac{\mathit{\text{E}}_i{}^{n+1}-\mathit{\text{E}}_i{}^n}{\mathrm{\Delta }t}={\epsilon }^{-1}\mathbf{M}_i{}^{-1}\left({\mathbf{S}}_i\mathit{\text{H}}_i{}^{n+\frac{1}{2}}+f_e{}_i(\mathit{\text{u}})\right)\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill =\mathbf{A}_e{}_i\mathit{\text{u}}(t),\hfill & & & & & & & & & & & & & & & & & & & \end{array}$$

(7)

where $\mathrm{\Delta }t$ is the size of the time step, and $f\left(\mathit{\text{u}}\right)$ denotes the flux term:

$$\begin{array}{cccccccccccccccccccc}\hfill f_h(\mathit{\text{u}})=-\sum _{p=1}^4\left(k_e\mathbf{F}{}_{{\mathrm{\Gamma }}_p}{}_i(\mathit{\text{E}}_j{}^n-\mathit{\text{E}}_i{}^n)\right)+\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \sum _{p=1}^4\left(v_h\mathbf{G}{}_{{\mathrm{\Gamma }}_p}{}_i(\mathit{\text{H}}_j{}^{n-\frac{1}{2}}-\mathit{\text{H}}_i{}^{n-\frac{1}{2}})\right),\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill f_e(\mathit{\text{u}})=\sum _{p=1}^4\left(k_h\mathbf{F}{}_{{\mathrm{\Gamma }}_p}{}_i(\mathit{\text{H}}_j{}^{n+\frac{1}{2}}-\mathit{\text{H}}_i{}^{n+\frac{1}{2}})\right)-\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \sum _{p=1}^4\left(v_e\mathbf{G}{}_{{\mathrm{\Gamma }}_p}{}_i(\mathit{\text{E}}_j{}^n-\mathit{\text{E}}_i{}^n)\right).\hfill & & & & & & & & & & & & & & & & & & & \end{array}$$

(8)

The half step between markedasmathbold italicE and markedasmathbold italicH proves that (8) has the second-order time accuracy [8]. The CFL condition constrains the $\mathrm{\Delta }t$ size [8, 9]:

$$\begin{array}{cccccccccccccccccccc}\hfill 0<\frac{\mathrm{\Delta }{t}^2}{4}max(\lambda \left(\mathbf{A}_e{}_i\right),\lambda \left(\mathbf{A}_h{}_i\right))<1,\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathrm{\Delta }{t}_{LF}\le \frac{2}{\sqrt{max(\lambda \left(\mathbf{A}_e{}_i\right),\lambda \left(\mathbf{A}_h{}_i\right))}},i\in \mathrm{\Omega }\hfill & & & & & & & & & & & & & & & & & & & \end{array}.$$

(9)

where $\lambda$ is the eigenvalue, and $\mathrm{\Delta }{t}_{LF}$ is the time step size of LF scheme.

B. Alternate boundary LTS scheme

The LF-LTS scheme can be used to expand (8) into the following:

$$\{\begin{array}{cccccccccccccccccccc}\hfill \frac{\mathit{\text{H}}_i{}^{n+\frac{1}{2}}-\mathit{\text{H}}_i{}^{n-\frac{1}{2}}}{\mathrm{\Delta }t}=\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathbf{A}_h{}_i{[\mathit{\text{E}}_i{}^n,\mathit{\text{E}}_j{}^n,\mathit{\text{H}}_i{}^{n-\frac{1}{2}},\mathit{\text{H}}_j{}^{n-\frac{1}{2}}]}^T\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \frac{\mathit{\text{E}}_i{}^{n+1}-\mathit{\text{E}}_i{}^n}{\mathrm{\Delta }t}=\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathbf{A}_e{}_i{[\mathit{\text{E}}_i{}^n,\mathit{\text{E}}_j{}^n,\mathit{\text{H}}_i{}^{n+\frac{1}{2}},\mathit{\text{H}}_j{}^{n+\frac{1}{2}}]}^T\hfill & & & & & & & & & & & & & & & & & & & \end{array},i\in e_c$$

(10)

$$\{\begin{array}{cccccccccccccccccccc}\hfill \frac{\mathit{\text{H}}_i{}^{n+\frac{1}{6}}-\mathit{\text{H}}_i{}^{n-\frac{1}{6}}}{\mathrm{\Delta }t\mathrm{/}3}=\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathbf{A}_h{}_i{[\mathit{\text{E}}_i{}^n,\mathit{\text{E}}_j{}^n,\mathit{\text{H}}_i{}^{n-\frac{1}{6}},\mathit{\text{H}}_j{}^{n-\frac{1}{6}}]}^T\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \frac{\mathit{\text{E}}_i{}^{n+\frac{1}{3}}-\mathit{\text{E}}_i{}^n}{\mathrm{\Delta }t\mathrm{/}3}=\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathbf{A}_h{}_i{[\mathit{\text{E}}_i{}^n,\mathit{\text{E}}_j{}^n,\mathit{\text{H}}_i{}^{n+\frac{1}{6}},\mathit{\text{H}}_j{}^{n+\frac{1}{6}}]}^T\hfill & & & & & & & & & & & & & & & & & & & \end{array},i\in e_f,$$

(11)

where $e_f$ and $e_c$ are the fine and coarse regions, respectively.

When the iteration steps of the fine region become an integer multiple of 3, the coarse region updates fields. According to the difference in the element scale, the computing domain can be classified into multiple LTS levels. The element of each level is coarse and fine for the upper and lower levels, respectively. Elements in the same LTS level ${\mathrm{\Omega }}_{LTS}$ have the same time step $\mathrm{\Delta }{t}_{LTS}$. A report [15] proved that

$$\begin{array}{cccccccccccccccccccc}\hfill 0<\frac{\min {\left(\mathrm{\Delta }{t}_{LTS}\right)}^2}{4}max\left(\lambda \left({\mathbf{A}}_p\right)\right)<1,\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mathbf{A}}_p={\mathbf{A}}_{\mathrm{\Omega }}-\frac{\min {\left(\mathrm{\Delta }{t}_{LTS}\right)}^2}{8}{\mathbf{A}}_{\mathrm{\Omega }}{\mathbf{A}}_{fine},\hfill & & & & & & & & & & & & & & & & & & & \end{array}$$

(12)

where ${\mathbf{A}}_{\mathrm{\Omega }}$ is the matrix established using elements in $\mathrm{\Omega }$, ${\mathbf{A}}_{fine}$ is the matrix established by employing elements in the fine region.

If markedasmathitalici is in the fine region and markedasmathitalicj is in the coarse region, (11) becomes

$$\{\begin{array}{cccccccccccccccccccc}\hfill \frac{\mathit{\text{H}}_i{}^{n+\frac{1}{6}}-\mathit{\text{H}}_i{}^{n-\frac{1}{6}}}{\mathrm{\Delta }t\mathrm{/}3}=\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mathbf{A}}_h{[\mathit{\text{E}}_i{}^n,\mathit{\text{E}}_j{}^n,\mathit{\text{H}}_i{}^{n-\frac{1}{2}},\mathit{\text{H}}_j{}^{n-\frac{1}{2}}]}^T\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \frac{\mathit{\text{E}}_i{}^{n+\frac{1}{3}}-\mathit{\text{E}}_i{}^n}{\mathrm{\Delta }t\mathrm{/}3}=\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mathbf{A}}_e{[\mathit{\text{E}}_i{}^n,\mathit{\text{E}}_j{}^n,\mathit{\text{H}}_i{}^{n-\frac{1}{2}},\mathit{\text{H}}_j{}^{n-\frac{1}{2}}]}^T\hfill & & & & & & & & & & & & & & & & & & & \end{array},i\in e_f,$$

(13)

The field obtained using j is not at the exact time step and reduces the accuracy and stability; thus, the time accuracy of (13) has the first order. In general, efficient methods, such as introduction of interpolation or exact iteration in an interface, are used to refine this scheme [6, 9–11], which increases the computing time and requires additional memory.

To avoid introducing interpolation, the boundary of fine and coarse elements is reconstructed so that the difference scale of coarse elements on the boundary is changed to be consistent with that of fine elements. The changed elements are alternate boundary (AB) elements, as displayed as Fig. 1.

Figure 1: Common case of alternate boundary method in which the yellow elements are alternate boundary elements that belong to the coarse region.

Considering (5), a different scheme of AB elements can be transformed into the following expression as the second-order scheme

$$\begin{array}{cccccccccccccccccccc}\hfill {\mathit{\text{u}}}^{n+\frac{1}{3}}\simeq \mathit{\text{u}}\left(t+\frac{\mathrm{\Delta }t}{3}\right)=\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathit{\text{u}}\left(t\right)+\frac{\mathrm{\Delta }t}{3}\frac{\partial \mathit{\text{u}}\left(t\right)}{\partial t}+\frac{1}{2}\left({\frac{\mathrm{\Delta }t}{3}}^2\right)\frac{{\partial }^2\mathit{\text{u}}\left(t\right)}{\partial t^2}+O(\mathrm{\Delta }{t}^3).\hfill & & & & & & & & & & & & & & & & & & & \end{array}$$

(14)

(14) provides the exact field of the coarse region and retains the second-order time accuracy in fine regions.

Different from the 2-D nodal basis function reported in [8, 9], the contribution of 3-D hierarchical vector basis functions to the upwind flux is mainly concentrated in the adjacent face and edges, and the orders are 0.5, 1.5, and 2, which render the spatial accuracy consistent with the time accuracy. The reciprocity of linear equations with the same order in space and time was confirmed in [16]. (14) of element i in the coarse region becomes:

$$\frac{\mathit{\text{u}}_i{}^{n+\frac{1}{3}}-\mathit{\text{u}}_i{}^n}{\mathrm{\Delta }t\mathrm{/}3}=\mathbf{A}\mathit{\text{u}}_i{}^n+\frac{1}{2}\mathbf{A}\left[\frac{\mathrm{\Delta }t}{2}\mathbf{A}\mathit{\text{u}}_i{}^{n-\frac{1}{2}}\right]+O(\mathrm{\Delta }{t}^3),$$

(15)

rewriting the equation in the LF scheme, we obtain

$$\begin{array}{cccccccccccccccccccc}\hfill \frac{\mathit{\text{H}}_i{}^{n+\frac{1}{6}}-\mathit{\text{H}}_i{}^{n-\frac{1}{6}}}{\mathrm{\Delta }t\mathrm{/}3}={\mathbf{A}}_h{[\mathit{\text{E}}_i{}^n,\mathit{\text{H}}_i{}^{n-\frac{1}{6}},\mathit{\text{E}}_j{}^n,\mathit{\text{H}}_j{}^{n-\frac{1}{6}}]}^T\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill +\frac{\mathrm{\Delta }t}{4}{\mathbf{A}}_h{\mathbf{A}}_e{[\mathit{\text{E}}_i{}^n,\mathit{\text{H}}_i{}^{n-\frac{1}{6}},\mathit{\text{E}}_j{}^n,\mathit{\text{H}}_j{}^{n-\frac{1}{6}}]}^T+O(\mathrm{\Delta }{t}^3).\hfill & & & & & & & & & & & & & & & & & & & \end{array}$$

(16)

(16) indicates that only two iterations of the AB element are required to obtain a field with high-order accuracy, according to (15). (16) also maintains the characteristics of the explicit scheme; hence, the AB elements must meet the CFL conditions of the fine region.

C. Computation work of ABLTS

It’s obvious that (16) is able to suppress errors in the coarse elements, and there will be an exact time step from the adjacent fine elements when (16) is applied in the upwind flux. Therefore, AB elements only need to maintain the same update scheme as adjacent fine elements to ensure that the accuracy of adjacent coarse elements is maintained at the second-order.

This feature makes the computation work of ABLTS method simpler than interpolation methods in [15]. For example, when the problem has 3 classes, the computation work of ABLTS is shown in Fig. 2.

The update for electric field of AB elements at the exact time is similar to (16):

$$\begin{array}{cccccccccccccccccccc}\hfill \frac{\mathit{\text{E}}_i{}^{n+\frac{1}{3}}-\mathit{\text{E}}_i{}^n}{\mathrm{\Delta }t\mathrm{/}3}={\mathbf{A}}_e{[\mathit{\text{E}}_i{}^n,\mathit{\text{H}}_i{}^{n+\frac{1}{6}},\mathit{\text{E}}_j{}^n,\mathit{\text{H}}_j{}^{n+\frac{1}{6}}]}^T\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill +\frac{\mathrm{\Delta }t}{4}{\mathbf{A}}_h{\mathbf{A}}_e{[\mathit{\text{E}}_i{}^n,\mathit{\text{H}}_i{}^{n+\frac{1}{6}},\mathit{\text{E}}_j{}^n,\mathit{\text{H}}_j{}^{n+\frac{1}{6}}]}^T+O(\mathrm{\Delta }{t}^3).\hfill & & & & & & & & & & & & & & & & & & & \end{array}$$

(17)

Figure 2: Computation work of ABLTS in the LF scheme, The time step size has $\mathrm{\Delta }{t}_2=3\mathrm{\Delta }{t}_1$; the red dotted lines are fields provided by Eq. (16).

From (16) and (17), the stability of the ABLTS method remains explicitly dependent on $\mathrm{\Delta }t$, and the upwind flux ensures equation convergence. The order of time in AB elements is the same as second-order interpolation. Therefore, the time step size of ABLTS follows (12).

III. NUMERICAL RESULTS AND DISCUSSION

In this section, the simulation of a resonator cavity and multi-layer antenna is presented to show the stability and efficiency of the ABLTS method.

A. Resonator cavity

First, the characteristics of a resonator cavity are analyzed to confirm the stability and accuracy of the proposed method. The cavity is 1 m × 1 m × 1 m and is terminated using a perfect electronic conductor (PEC) boundary. The interior of the resonant cavity is filled with air, so that the resonant frequency obtained using the analytical solution is 212.132 MHz. The cavity is partitioned to A, B, and C regions, and the mesh size ratio of these regions is set to 1:3:9 for showing the multi-scale situation.

Figure 3: Geometry and meshes of the resonator cavity.

To verify the stability and accuracy of the ABLTS method, three control groups are used for testing. The mesh size of region A in these control groups were 10, 15, and 20 mm. After mesh generation, the numbers of elements of these groups are 75,351, 35,679, and 17,114 tetrahedrons; the number of degree of freedoms (Dofs) are 3,014,040, 1,427,160, 684,560. A point source is excited at the center of the cavity with modulated Gaussian pulse from 100 to 300 MHz. The DGTD method with LF, LF-LTS presented in [9], and ABLTS methods are used to simulate the wave propagation of 1000 ns in the cavity. The minimum time step size of the above methods is 37.562 ps. Fig. 3 shows the geometry of the cavity.

Table 1 presents the performance comparison of the relative error and solution time of LF, LTS, and ABLTS methods, and the relative error is compared with the result of analytical solution of 212.132 MHz (absolute error/analytical solution). For the same mesh and excitation and solution times, the accuracy of the ABLTS method is the same as that of the LF method and higher than that of the LF-LTS method. The speed of the ABLTS method is similar to that of the LF-LTS method and considerably higher than that of the LF method. In this example, the spatial discretization at a maximum scale of 0.127 wavelengths and following eq. (9) and eq. (12) to estimate the time step, three groups of experiments at different scales all maintained stable working for a long time as Fig. 4 shows, demonstrating the stability of the proposed method. The electric field of the cavity is shown in Fig. 5.

Table 1: Performance comparison among different methods

Group	Performance	LF	LF-LTS	ABLTS
Group 1	Time (min)	315.4	90.64	90.75
	Speedup	$\backslash$	3.48	3.48
	Freq (MHz)	212.106	212.011	212.105
	Relative error	0.01%	0.06%	0.01%
Group 2	Time (min)	69.12	22.35	22.07
	Speedup	$\backslash$	3.09	3.08
	Freq (MHz)	212.095	212.003	212.093
	Relative error	0.02%	0.06%	0.02%
Group 3	Time (min)	27.51	9.09	9.1
	Speedup	$\backslash$	3.02	3.02
	Freq (MHz)	212.039	211.991	212.039
	Relative error	0.04%	0.07%	0.04%

Figure 4: Nearfield of the proposed method.

Figure 5: Electric field (V/m) of the cavity at 28 ns, 40 ns, 44 ns, and 54 ns.

B. Multi-layer microstrip antenna array

A multi-layer microstrip antenna array is analyzed using the DGTD-ABLTS method to show its efficiency and accuracy. The array had four layers (Fig. 6); the first layer of the upper part is made of Rogers RO3203, whose permittivity is 3.02. The second and forth layers are composed of Arlon CuClad 250GT, whose permittivity is 2.5. The third layer is a perfect electric conductor. The computation domain is terminated by the Silver-Muller absorb boundary condition. Four lump ports excited the feed network through a modulated Gaussian pulse from 13.5 GHz to 17.5 GHz. The number of tetrahedrons is 1,735,017, the number of Dofs is 69,400,680, and the ABLTS method decomposes meshes into three levels. The numbers of tetrahedron of each LTS levels are 42,264, 408,453, and 1,731,841. The pulse propagated in the model in 4 ns.

Figure 6: Geometry of the multi-layer microstrip antenna array.

Figure 7: Electric field (V/m) in the XOZ plane at 0.8 ns.

Figure 8: GainTotal in YOZ plane calculated by the proposed method and HFSS.

Table 2: Performance of the DGTD with ABLTS

	Peak	Solution	$\mathit{\delta }$	Speedup
	Memory	Time	(Compare	(Compare
	(MB)	(min)	with LF)	with LF)
LF	29,264.1	416.3	$\backslash$	$\backslash$
LTS	29,266.9	54.2	0.0937	7.68
ABLTS	29,267.3	54.7	0.0815	7.61

Figure 7 shows the electric field distribution in the XOZ plane at 0.8 ns, and Fig. 6 presents the Gain Total in the YOZ plane in comparison with the result of HFSS. Furthermore, Table 2 presents the computational performance of the LF-DGTD and proposed methods, in which $\delta$ denotes the relative error. All the methods are run on CPUs with 240 processors. The computation efficiency of the DGTD method considerably improved with ABLTS and the memory and accuracy are almost unchanged.

IV. CONCLUSIONS

An ABLTS method of DGTD is introduced for electromagnetic simulation of antennas. The DGTD-ABLTS method is a noninterpolation local time stepping scheme based on the leapfrog integration scheme, which reduces the error of the LTS method by two iterations with the second-order accuracy. The simulation of a resonator cavity and multi-layer antenna array proved the accuracy and efficiency of the DGTD-ABLTS method. The DGTD-ABLTS method is efficient for the large-scale electromagnetic simulation of antennas.

ACKNOWLEDGMENT

This work was supported by the Key Research and Development Program of Shaanxi(2023-ZDLGY-09, 2021GXLH-02), and the Fundamental Research Funds for the Central Universities (QTZX22160), and the Special Funds of CEPREI (23Z013).

REFERENCES

[1] J. S. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, New York, NY, USA: Springer, 2008.

[2] A. Bossavit, “Whitney forms: a class of finite elements for three-dimensional computations in electromagnetism,” IEE Proceedings A Physical Science, Measurement and Instrumentation, Management and Education, Reviews, vol. 135, no. 8, pp. 493-500, 1988.

[3] X. Zhao, Z. Lin, Y. Zhang, S. Ting, and T. K. Sarkar, “Parallel Hybrid method of HOMoM–MLFMA for analysis of large antenna arrays on an electrically large platform,” IEEE Transactions on Antennas and Propagation, vol. 64, no. 12, pp. 5501-5506, Dec. 2016, doi:10.1109/TAP.2016.2621029.

[4] L. Zhao, G. Chen, W. Yu, and J. M. Jin, “A fast waveguide port parameter extraction technique for the DGTD method,” IEEE Antennas and Wireless Propagation Letters, vol. 16, pp. 2659-2662, 2017, doi:10.1109/LAWP.2017.2740298.

[5] S. Dosopoulos and J. F. Lee, “Interior penalty discontinuous Galerkin finite element method for the time-dependent first order Maxwell’s equations,” IEEE Transactions on Antennas and Propagation, vol. 58, no. 12, pp. 4085-4090, Dec. 2010, doi:10.1016/j.jcp.2010.07.036.

[6] J. Chen and Q. Liu, “Discontinuous Galerkin time-domain methods for multiscale electromagnetic simulations: A review,” Proceedings of the IEEE, vol. 101, no. 2, pp. 242-254, Feb. 2013, doi: 10.1109/JPROC.2012.2219031.

[7] L. Zhao, G. Chen, W. Yu, and J. M. Jin, “A fast waveguide port parameter extraction technique for the DGTD method,” IEEE Antennas and Wireless Propagation Letters, vol. 16, pp. 2659-2662, 2017, doi:10.1109/LAWP.2017.2740298.

[8] M. J. Grote and M. Teodora. “Explicit local time-stepping methods for time-dependent wave propagation,” arXiv: Numerical Analysi, 2012, doi:10.48550/arXiv.1205.0654.

[9] E. Montseny, S. Pernet, X. Ferrières, and G. Cohen, “Dissipative terms and local time-stepping improvements in a spatial high order Discontinuous Galerkin scheme for the time-domain Maxwell’s equations,” Journal of Computational Physics, vol. 227, no. 14, pp. 6795-6820, 2008.

[10] L. Fezoui, S. Lanteri, S. Lohrengel, and S. Piperno, “Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes,” ESAIM, Mathematical Modelling and Numerical Analysis, vol. 39, no. 6, pp. 1149-1176, 2005.

[11] J. Alvarez, L. D. Angulo, A. R. Bretones, C. M. Coevorden and S. G. Garcia, “Efficient antenna modeling by DGTD: Leap-frog discontinuous Galerkin timedomain method,” IEEE Antennas and Propagation Magazine, vol. 57, no. 3, pp. 95-106, June 2015, doi:10.1109/MAP.2015.2437279.

[12] M. Li, Q. Wu, Z. Lin, Y. Zhang, and X. Zhao, “Novel parallelization of discontinuous galerkin method for transient electromagnetics simulation based on sunway supercomputers,” Applied Computational Electromagnetics Society (ACES) Journal, vol. 37, no. 7, pp. 795-804, Sep. 2021, doi:10.13052/2022.ACES.J.370706.

[13] M. Li, Q. Wu, Z. Lin, Y. Zhang, and X. Zhao, “A minimal round-trip strategy based on graph matching for parallel dgtd method with local Time-stepping,” IEEE Antennas and Wireless Propagation Letters, vol. 22, no. 2, pp. 243-247, Feb. 2023, doi:10.1109/LAWP.2022.3208010.

[14] Z. G. Ban, Y. Shi, and P. Wang, “Advanced parallelism of DGTD method with local time stepping based on novel MPI + MPI unified parallel algorithm,” IEEE Transactions on Antennas and Propagation, vol. 70, no. 5, pp. 3916-3921, 2022, doi:10.1109/TAP.2021.3137455.

[15] G. Cohen, X. Ferrieres, and S. Pernet, “A spatial high-order hexahedral discontinuous Galerkin method to solve Maxwell’s equations in time domain,” Journal of Computational Physics, vol. 217, no. 2, pp. 340-363, 2006, doi:10.1016/j.jcp.2006.01.004.

[16] S. Zuo, Z. Lin, D. G. Doñoro, Y. Zhang, and X. Zhao, “A parallel direct domain decomposition solver based on schur complement for electromagnetic finite element analysis,” IEEE Antennas and Wireless Propagation Letters, vol. 20, no. 4, pp. 458-462, Apr. 2021, doi:10.1109/LAWP.2021.3053566.

BIOGRAPHIES

Minxuan Li was born in Changzhi, Shanxi, China, in 1994. He received the B.S. degree in electronic information engineering from Harbin Institute of Technology in 2016. He is currently pursuing the Ph.D. degree with Xidian University. He joined the China Institute of Electronic Product Reliability and Environmental Testing in 2023. His research interests include channel modeling, parallel computation, and transient electromagnetic analysis.

Qingkai Wu was born in Yangzhou, Zhejiang, China, in 1997. He received the B.S. degree in electronic science and technology from Xidian University in 2020. He is currently pursuing the Ph.D. degree with Xidian University. His current research interests include transient electromagnetic analysis.

Zhongchao Lin was born in Hubei, China, in 1988. He received the B.S. and Ph.D. degrees from Xidian University, in 2011 and 2016, respectively. He joined Xidian University, in 2016, as a post doctoral fellow, where he was lately promoted as an associate professor. His research interests include large-scale computational electromagnetic, scattering, and radiation electromagneticanalysis.

YU ZHANG received the B.S., M.S., and Ph.D. degrees from Xidian University, in 1999, 2002, and 2004, respectively. In 2004, he joined Xidian University as a faculty member. He was a visiting scholar and an adjunct professor with Syracuse University, from 2006 to 2009. As a principal investigator, he works on projects, including the project of NSFC. He has authored four books: Parallel Computation in Electromagnetics (Xidian University Press, 2006), Parallel Solution of Integral Equation-Based EM Problems in the Frequency Domain (Wiley IEEE, 2009), Time and Frequency Domain Solutions of EM Problems Using Integral Equations and a Hybrid Methodology (Wiley, 2010), and Higher Order Basis Based Integral Equation Solver (Wiley, 2012), as well as more than 100 journal articles and 40 conferencepapers.

Xunwang Zhao was born in Shanxi, China, in 1983. He received the B.S. and Ph.D. degrees from Xidian University, in 2004 and 2008, respectively. He joined Xidian University, in 2008, as a faculty member, where he was lately promoted as a full professor. He was a visiting scholar with Syracuse University, from December 2008 to April 2009. As a principal investigator, he works on several projects, including the project of NSFC. His research interests include computational electromagnetic and electromagnetic scattering analysis.

ACES JOURNAL, Vol. 38, No. 8, 587–593
doi: 10.13052/2023.ACES.J.380805
© 2023 River Publishers

ABSTRACT

I. INTRODUCTION

II. PRINCIPLE AND FORMULATIONS

A. Numerical discretization of LTS-DGTD

B. Alternate boundary LTS scheme

C. Computation work of ABLTS

III. NUMERICAL RESULTS AND DISCUSSION

A. Resonator cavity

B. Multi-layer microstrip antenna array

IV. CONCLUSIONS

ACKNOWLEDGMENT

REFERENCES

BIOGRAPHIES