A Dissipation Theory for Three-Dimensional FDTD with Application to Stability Analysis and Subgridding

TL;DR

This paper introduces a dissipation-based framework for 3-D FDTD that ensures stability in complex setups, enabling modular design of stable schemes like a new subgridding method with low reflections and flexible refinement.

Contribution

It presents a novel dissipation theory for 3-D FDTD that guarantees stability by design, simplifying the creation of stable, modular electromagnetic simulation schemes.

Findings

Proposes a dissipation framework inspired by energy conservation.

Derives a stable subgridding scheme with low reflections.

Enables modular stability enforcement for complex FDTD setups.

Abstract

The finite-difference time-domain (FDTD) algorithm is a popular numerical method for solving electromagnetic problems. FDTD simulations can suffer from instability due to the explicit nature of the method. Stability enforcement can be particularly challenging in scenarios where a setup is composed of multiple components, such as grids of different resolution, advanced boundary conditions, reduced-order models, and lumped elements. We propose a dissipation theory for 3-D FDTD inspired by the principle of energy conservation. We view the FDTD update equations for a 3-D region as a dynamical system, and show under which conditions the system is dissipative. By requiring each component of an FDTD-like scheme to be dissipative, the stability of the overall coupled scheme follows by construction. The proposed framework enables the creation of provably stable schemes in an easy and modular fashion, since conditions are imposed on the individual components, rather than on the overall coupled scheme as in existing approaches. With the proposed framework, we derive a new subgridding scheme with guaranteed stability, low reflections, support for material traverse and arbitrary (odd) grid refinement ratio.

Tables

Table 1. Table I: Computational Cost for Different Discretization in the Meta-Screen Example of Sec.  VI-C .

Method	Number of primary cells	Simulation time	Speed-up
Coarse	00,281,232	00.18 hours	524.8
Fine	41,341,104	95.81 hoursa	reference
Nonuniform	04,065,600	11.46 hours	8.4
Proposed	01,323,672	05.80 hours	16.5

Equations

$$\Delta x\Delta y\Delta z\left(\frac{\varepsilon_{x}}{\Delta t}\!+\!\frac{\sigma_{x}}{2}\right)E_{x}\big{|}_{i+\frac{1}{2},j,k}^{n+1}\\ =\Delta x\Delta y\Delta z\left(\frac{\varepsilon_{x}}{\Delta t}\!-\!\frac{\sigma_{x}}{2}\right)E_{x}\big{|}_{i+\frac{1}{2},j,k}^{n}\\ +\Delta x\Delta zH_{z}\big{|}_{i+\frac{1}{2},j+\frac{1}{2},k}^{n+\frac{1}{2}}-\Delta x\Delta yH_{y}\big{|}_{i+\frac{1}{2},j,k+\frac{1}{2}}^{n+\frac{1}{2}}\\ -\Delta x\Delta zH_{z}\big{|}_{i+\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2}}+\Delta x\Delta yH_{y}\big{|}_{i+\frac{1}{2},j,k-\frac{1}{2}}^{n+\frac{1}{2}}\,,$$

$$\Delta x\frac{\Delta y}{2}\Delta z\left(\frac{\varepsilon_{x}}{\Delta t}\!+\!\frac{\sigma_{x}}{2}\right)E_{x}\big{|}_{i+\frac{1}{2},1,k}^{n+1}\\ =\Delta x\frac{\Delta y}{2}\Delta z\left(\frac{\varepsilon_{x}}{\Delta t}\!-\!\frac{\sigma_{x}}{2}\right)E_{x}\big{|}_{i+\frac{1}{2},1,k}^{n}\\ +\Delta x\Delta zH_{z}\big{|}_{i+\frac{1}{2},\frac{3}{2},k}^{n+\frac{1}{2}}-\Delta x\frac{\Delta y}{2}H_{y}\big{|}_{i+\frac{1}{2},1,k+\frac{1}{2}}^{n+\frac{1}{2}}\\ -\Delta x\Delta zU_{z}\big{|}_{i+\frac{1}{2},1,k}^{n+\frac{1}{2}}+\Delta x\frac{\Delta y}{2}H_{y}\big{|}_{i+\frac{1}{2},1,k-\frac{1}{2}}^{n+\frac{1}{2}}\,.$$

$$\Delta x\frac{\Delta y}{2}\frac{\Delta z}{2}\left(\frac{\varepsilon_{x}}{\Delta t}\!+\!\frac{\sigma_{x}}{2}\right)E_{x}\big{|}_{i+\frac{1}{2},1,1}^{n+1}\\ =\Delta x\frac{\Delta y}{2}\frac{\Delta z}{2}\left(\frac{\varepsilon_{x}}{\Delta t}\!-\!\frac{\sigma_{x}}{2}\right)E_{x}\big{|}_{i+\frac{1}{2},1,1}^{n}\\ +\Delta x\frac{\Delta z}{2}H_{z}\big{|}_{i+\frac{1}{2},\frac{3}{2},1}^{n+\frac{1}{2}}-\Delta x\frac{\Delta y}{2}H_{y}\big{|}_{i+\frac{1}{2},1,\frac{3}{2}}^{n+\frac{1}{2}}\\ -\Delta x\frac{\Delta z}{2}U_{z}\big{|}_{i+\frac{1}{2},1,1}^{n+\frac{1}{2}}+\Delta x\frac{\Delta y}{2}U_{y}\big{|}_{i+\frac{1}{2},1,1}^{n+\frac{1}{2}}\,.$$

$$\Delta x\Delta y\Delta z\frac{\mu_{x}}{\Delta t}H_{x}\big{|}_{i,j+\frac{1}{2},k+\frac{1}{2}}^{n+\frac{1}{2}}=\Delta x\Delta y\Delta z\frac{\mu_{x}}{\Delta t}H_{x}\big{|}_{i,j+\frac{1}{2},k+\frac{1}{2}}^{n-\frac{1}{2}}\\ -\Delta x\Delta zE_{z}\big{|}_{i,j+1,k+\frac{1}{2}}^{n}+\Delta x\Delta yE_{y}\big{|}_{i,j+\frac{1}{2},k+1}^{n}\\ +\Delta x\Delta zE_{z}\big{|}_{i,j,k+\frac{1}{2}}^{n}-\Delta x\Delta yE_{y}\big{|}_{i,j+\frac{1}{2},k}^{n}\,,$$

$$\frac{\Delta x}{2}\Delta y\Delta z\frac{\mu_{x}}{\Delta t}H_{x}\big{|}_{1,j+\frac{1}{2},k+\frac{1}{2}}^{n+\frac{1}{2}}=\frac{\Delta x}{2}\Delta y\Delta z\frac{\mu_{x}}{\Delta t}H_{x}\big{|}_{1,j+\frac{1}{2},k+\frac{1}{2}}^{n-\frac{1}{2}}\\ -\frac{\Delta x}{2}\Delta zE_{z}\big{|}_{1,j+1,k+\frac{1}{2}}^{n}+\frac{\Delta x}{2}\Delta yE_{y}\big{|}_{1,j+\frac{1}{2},k+1}^{n}\\ +\frac{\Delta x}{2}\Delta zE_{z}\big{|}_{1,j,k+\frac{1}{2}}^{n}-\frac{\Delta x}{2}\Delta yE_{y}\big{|}_{1,j+\frac{1}{2},k}^{n}\,,$$

$$\begin{gathered} {\mathbf{D}}_{l_x^{\prime }}{\mathbf{D}}_{A_x}\frac{{\mathbf{D}}_{{\mu }_x}}{\Delta t}{\mathbf{H}}_x^{n+\frac{1}{2}}={\mathbf{D}}_{l_x^{\prime }}{\mathbf{D}}_{A_x}\frac{{\mathbf{D}}_{{\mu }_x}}{\Delta t}{\mathbf{H}}_x^{n-\frac{1}{2}} \\ +{\mathbf{D}}_{l_x^{\prime }}{\mathbf{G}}_{zy}{\mathbf{D}}_{l_y}{\mathbf{E}}_y^n-{\mathbf{D}}_{l_x^{\prime }}{\mathbf{G}}_{yz}{\mathbf{D}}_{l_z}{\mathbf{E}}_z^n, \end{gathered}$$

$$\begin{gathered} {\mathbf{D}}_{l_y^{\prime }}{\mathbf{D}}_{A_y}\frac{{\mathbf{D}}_{{\mu }_y}}{\Delta t}{\mathbf{H}}_y^{n+\frac{1}{2}}={\mathbf{D}}_{l_y^{\prime }}{\mathbf{D}}_{A_y}\frac{{\mathbf{D}}_{{\mu }_y}}{\Delta t}{\mathbf{H}}_y^{n-\frac{1}{2}} \\ -{\mathbf{D}}_{l_y^{\prime }}{\mathbf{G}}_{zx}{\mathbf{D}}_{l_x}{\mathbf{E}}_x^n+{\mathbf{D}}_{l_y^{\prime }}{\mathbf{G}}_{xz}{\mathbf{D}}_{l_z}{\mathbf{E}}_z^n, \end{gathered}$$

$$\begin{gathered} {\mathbf{D}}_{l_z^{\prime }}{\mathbf{D}}_{A_z}\frac{{\mathbf{D}}_{{\mu }_z}}{\Delta t}{\mathbf{H}}_z^{n+\frac{1}{2}}={\mathbf{D}}_{l_z^{\prime }}{\mathbf{D}}_{A_z}\frac{{\mathbf{D}}_{{\mu }_z}}{\Delta t}{\mathbf{H}}_z^{n-\frac{1}{2}} \\ +{\mathbf{D}}_{l_z^{\prime }}{\mathbf{G}}_{yx}{\mathbf{D}}_{l_x}{\mathbf{E}}_x^n-{\mathbf{D}}_{l_z^{\prime }}{\mathbf{G}}_{xy}{\mathbf{D}}_{l_y}{\mathbf{E}}_y^n, \end{gathered}$$

$$\begin{gathered} {\mathbf{D}}_{l_x}{\mathbf{D}}_{A_x^{\prime }}\left(\frac{{\mathbf{D}}_{{\epsilon }_x}}{\Delta t}+\frac{{\mathbf{D}}_{{\sigma }_x}}{2}\right){\mathbf{E}}_x^{n+1}={\mathbf{D}}_{l_x}{\mathbf{D}}_{A_x^{\prime }}\left(\frac{{\mathbf{D}}_{{\epsilon }_x}}{\Delta t}-\frac{{\mathbf{D}}_{{\sigma }_x}}{2}\right){\mathbf{E}}_x^n \\ +{\mathbf{D}}_{l_x}{\mathbf{G}}_{zx}^T{\mathbf{D}}_{l_y^{\prime }}{\mathbf{H}}_y^{n+\frac{1}{2}}-{\mathbf{D}}_{l_x}{\mathbf{G}}_{yx}^T{\mathbf{D}}_{l_z^{\prime }}{\mathbf{H}}_z^{n+\frac{1}{2}} \\ +{\mathbf{D}}_{l_x}{\mathbf{B}}_{x,B}{\mathbf{D}}_{l_y^{\prime },B}{\mathbf{U}}_{y,B}^{n+\frac{1}{2}}+{\mathbf{D}}_{l_x}{\mathbf{B}}_{x,N}{\mathbf{D}}_{l_z^{\prime },N}{\mathbf{U}}_{z,N}^{n+\frac{1}{2}} \\ -{\mathbf{D}}_{l_x}{\mathbf{B}}_{x,T}{\mathbf{D}}_{l_y^{\prime },T}{\mathbf{U}}_{y,T}^{n+\frac{1}{2}}-{\mathbf{D}}_{l_x}{\mathbf{B}}_{x,S}{\mathbf{D}}_{l_z^{\prime },S}{\mathbf{U}}_{z,S}^{n+\frac{1}{2}}, \end{gathered}$$

$${\mathbf{D}}_{l^{\prime }}{\mathbf{D}}_A\frac{{\mathbf{D}}_{\mu }}{\Delta t}{\mathbf{H}}^{n+\frac{1}{2}}={\mathbf{D}}_{l^{\prime }}{\mathbf{D}}_A\frac{{\mathbf{D}}_{\mu }}{\Delta t}{\mathbf{H}}^{n-\frac{1}{2}}-{\mathbf{D}}_{l^{\prime }}\mathbf{C}{\mathbf{D}}_{\text{l}}{\mathbf{E}}^n,$$

$${\mathbf{E}}^n=\left[\begin{array}{c}{\mathbf{E}}_x^n\\ {\mathbf{E}}_y^n\\ {\mathbf{E}}_z^n\end{array}\right],{\mathbf{H}}^{n-\frac{1}{2}}=\left[\begin{array}{c}{\mathbf{H}}_x^{n-\frac{1}{2}}\\ {\mathbf{H}}_y^{n-\frac{1}{2}}\\ {\mathbf{H}}_z^{n-\frac{1}{2}}\end{array}\right].$$

$$\mathbf{C}=\left[\begin{array}{ccc}0& -{\mathbf{G}}_{zy}& {\mathbf{G}}_{yz}\\ {\mathbf{G}}_{zx}& 0& -{\mathbf{G}}_{xz}\\ -{\mathbf{G}}_{yx}& {\mathbf{G}}_{xy}& 0\end{array}\right]$$

$$\begin{gathered} {\mathbf{D}}_{\text{l}}{\mathbf{D}}_{A^{\prime }}\left(\frac{{\mathbf{D}}_{\epsilon }}{\Delta t}+\frac{{\mathbf{D}}_{\sigma }}{2}\right){\mathbf{E}}^{n+1}={\mathbf{D}}_{\text{l}}{\mathbf{D}}_{A^{\prime }}\left(\frac{{\mathbf{D}}_{\epsilon }}{\Delta t}-\frac{{\mathbf{D}}_{\sigma }}{2}\right){\mathbf{E}}^n \\ +{\mathbf{D}}_{\text{l}}{\mathbf{C}}^T{\mathbf{D}}_{l^{\prime }}{\mathbf{H}}^{n+\frac{1}{2}}+{\mathbf{D}}_{\text{l}}\mathbf{QP}{\mathbf{D}}_{l_U^{\prime }}{\mathbf{U}}^{n+\frac{1}{2}}, \end{gathered}$$

$$\mathbf{QP}=\left[\begin{array}{ccc}{\mathbf{B}}_x& 0& 0\\ 0& {\mathbf{B}}_y& 0\\ 0& 0& {\mathbf{B}}_z\end{array}\right],$$

$${\mathbf{B}}_x=\left[\begin{array}{cccc}{\mathbf{B}}_{x,B}& {\mathbf{B}}_{x,N}& -{\mathbf{B}}_{x,T}& -{\mathbf{B}}_{x,S}\end{array}\right],$$

$${\mathbf{B}}_y=\left[\begin{array}{cccc}{\mathbf{B}}_{y,W}& {\mathbf{B}}_{y,T}& -{\mathbf{B}}_{y,E}& -{\mathbf{B}}_{y,B}\end{array}\right],$$

$${\mathbf{B}}_z=\left[\begin{array}{cccc}{\mathbf{B}}_{z,S}& {\mathbf{B}}_{z,E}& -{\mathbf{B}}_{z,N}& -{\mathbf{B}}_{z,W}\end{array}\right].$$

$$(\mathbf{R}+\mathbf{F}){\mathbf{x}}^{n+1}$$

$${\mathbf{y}}^n$$

$${\mathbf{x}}^n=\left[\begin{array}{c}{\mathbf{E}}^n\\ {\mathbf{H}}^{n-\frac{1}{2}}\end{array}\right].$$

$${\mathbf{u}}^{n+\frac{1}{2}}=\left[\begin{array}{c}{\mathbf{U}}_{E_x}^{n+\frac{1}{2}}\\ {\mathbf{U}}_{E_y}^{n+\frac{1}{2}}\\ {\mathbf{U}}_{E_z}^{n+\frac{1}{2}}\end{array}\right],$$

$${\mathbf{y}}^n=\left[\begin{array}{c}{\mathbf{Y}}_{E_x}^n\\ {\mathbf{Y}}_{E_y}^n\\ {\mathbf{Y}}_{E_z}^n\end{array}\right],$$

$$\mathbf{R}=\left[\begin{array}{cc}{\mathbf{D}}_{\text{l}}{\mathbf{D}}_{A^{\prime }}\frac{{\mathbf{D}}_{\epsilon }}{\Delta t}& -\frac{1}{2}{\mathbf{D}}_{\text{l}}{\mathbf{C}}^T{\mathbf{D}}_{l^{\prime }}\\ -\frac{1}{2}{\mathbf{D}}_{l^{\prime }}\mathbf{C}{\mathbf{D}}_{\text{l}}& {\mathbf{D}}_{l^{\prime }}{\mathbf{D}}_A\frac{{\mathbf{D}}_{\mu }}{\Delta t}\end{array}\right],$$

$$\mathbf{F}=\left[\begin{array}{cc}{\mathbf{D}}_{\text{l}}{\mathbf{D}}_{A^{\prime }}\frac{{\mathbf{D}}_{\sigma }}{2}& -\frac{1}{2}{\mathbf{D}}_{\text{l}}{\mathbf{C}}^T{\mathbf{D}}_{l^{\prime }}\\ \frac{1}{2}{\mathbf{D}}_{l^{\prime }}\mathbf{C}{\mathbf{D}}_{\text{l}}& 0\end{array}\right],$$

$$\mathbf{B}=\left[\begin{array}{c}{\mathbf{D}}_{\text{l}}\mathbf{QP}{\mathbf{D}}_{l_U^{\prime }}\\ 0\end{array}\right],$$

$${\mathbf{L}}^T=\left[\begin{array}{cc}{\mathbf{Q}}^T& 0\end{array}\right].$$

$$\mathcal{E}\left({\mathbf{x}}^n\right)$$

$$\mathcal{E}\left(0\right)$$

$$\mathcal{E}\left({\mathbf{x}}^{n+1}\right)-\mathcal{E}\left({\mathbf{x}}^n\right)$$

$$\mathcal{E}\left({\mathbf{x}}^n\right)=\frac{\Delta t}{2}{\left({\mathbf{x}}^n\right)}^T\mathbf{R}{\mathbf{x}}^n,$$

$$s\left({\mathbf{y}}^n,{\mathbf{u}}^{n+\frac{1}{2}}\right)=\Delta t{\frac{\left({\mathbf{y}}^n+{\mathbf{y}}^{n+1}\right)}{2}}^T\mathbf{S}{\mathbf{u}}^{n+\frac{1}{2}}.$$

$$\mathbf{S}={\mathbf{D}}_{l_Y}{\mathbf{D}}_{l_U^{\prime }}\mathbf{P},$$

$$\begin{gathered} \mathcal{E}\left({\mathbf{x}}^n\right)=\frac{1}{2}({\mathbf{E}}^n{)}^T{\mathbf{D}}_{\text{l}}{\mathbf{D}}_{A^{\prime }}{\mathbf{D}}_{\epsilon }{\mathbf{E}}^n \\ +\frac{\Delta t}{2}{\left({\mathbf{H}}^{n-\frac{1}{2}}\right)}^T\left[{\mathbf{D}}_{l^{\prime }}{\mathbf{D}}_A\frac{{\mathbf{D}}_{\mu }}{\Delta t}{\mathbf{H}}^{n-\frac{1}{2}}-{\mathbf{D}}_{l^{\prime }}\mathbf{C}{\mathbf{D}}_{\text{l}}{\mathbf{E}}^n\right], \end{gathered}$$

$$\begin{gathered} \mathcal{E}\left({\mathbf{x}}^n\right)=\frac{1}{2}({\mathbf{E}}^n{)}^T{\mathbf{D}}_{\text{l}}{\mathbf{D}}_{A^{\prime }}{\mathbf{D}}_{\epsilon }{\mathbf{E}}^n \\ +\frac{1}{2}{\left({\mathbf{H}}^{n-\frac{1}{2}}\right)}^T{\mathbf{D}}_{l^{\prime }}{\mathbf{D}}_A{\mathbf{D}}_{\mu }{\mathbf{H}}^{n+\frac{1}{2}}. \end{gathered}$$

Full text

I Introduction

Stability is a major challenge in explicit numerical schemes for solving differential equations. The conventional finite-difference time-domain (FDTD) algorithm for solving Maxwell’s equations [1] is stable when the iteration time step is below the Courant-Friedrichs-Lewy (CFL) limit [2]. However, stability enforcement becomes much more difficult in the case of more advanced FDTD schemes such as, for example, those including locally-refined grids [3, 4, 5], reduced-order models [6, 7], and hybridizations of FDTD with integral equation methods [8] or circuit simulators [9]. Many of these advanced schemes can be seen as the coupling of multiple blocks, such as FDTD grids of different resolution, reduced-order models, lumped components, boundary conditions, and so on. While the stability of each individual block may be well understood, analyzing and enforcing the stability of the overall coupled scheme can be a daunting task.

Several techniques have been proposed for FDTD stability analysis, such as von Neumann analysis [10], the iteration method [2] and the energy method [11]. Unfortunately, von Neumann analysis cannot be applied in presence of inhomogeneous materials. The iteration method [2, 12] analyzes the eigenvalues of a matrix describing FDTD iterations. Even though an FDTD setup can often be divided into several subsystems, such as grids of different resolution, reduced-order models, interpolation conditions, and other blocks, the iteration matrix needs to be derived for the entire scheme, which can lead to lengthy derivations and make stability analysis very tedious. A similar issue arises in the case of the energy method [11], where one must write an expression for the stored energy in the entire domain, and show that the coupled scheme satisfies the principle of energy conservation.

In this paper, we present a dissipation theory for modular stability enforcement of complex 3-D FDTD systems, generalizing previous work in two dimensions [13]. The method is based on the theory of dissipative dynamical systems [14]. First, we write the FDTD update equations for an arbitrary inhomogeneous region in the form of a discrete-time dynamical system. The tangential magnetic and electric field are seen as input and output, respectively. Suitable discrete expressions for the energy stored inside the region, and for the energy absorbed through the boundaries are proposed and used to investigate under which conditions the system is dissipative. This result becomes the basis for a powerful framework to create new FDTD algorithms with guaranteed stability. By imposing each subsystem to to be dissipative, the overall coupled scheme formed by these subsystems is dissipative by construction, and thus stable. The main virtue of this approach is that the stability of the overall scheme follows automatically from conditions imposed on each subsystem individually. This makes stability analysis simpler and more modular, which is a remarkable advantage over existing approaches that require the analysis of the overall coupled scheme. This modularity also facilitates the generation of new schemes, since proving a given subsystem (e.g. a sophisticated boundary condition) to be dissipative allows coupling it to any other dissipative block, with no need for additional analysis. Finally, we use the theory to derive a stable subgridding scheme for 3-D FDTD that works for any odd refinement ratio and naturally supports traversal of material boundaries by the subgridding interface.

This paper is structured as follows. In Sec. II, we show how the update equations for an FDTD region can be cast into the form of a dynamical system with suitable inputs and outputs. In Sec. III, we provide dissipativity conditions for the system and show their relation to the CFL limit. In Sec. IV, we describe the modular method for enforcing FDTD stability using the proposed theory, and in Sec. V we apply this method to derive a stable subgridding scheme. Numerical examples are given in Sec. VI.

II Dynamical System Formulation of a 3-D FDTD Region

Consider the FDTD region shown in Fig. 1, which contains $N_{x}$, $N_{y}$, and $N_{z}$ primary cells in the $x$, $y$, and $z$ directions, respectively. Electric field is sampled on the primary grid edges at the integer time points, $n$, and magnetic field is sampled on the primary grid faces at the time points shifted by half of a time step, $n\!-\!(1/2)$ [2]. For simplicity, we assume that the region is discretized uniformly with primary cells of dimensions $\Delta x\!\times\!\Delta y\!\times\!\Delta z$. As shown in Fig. 1, we use cardinal directions in the $xy$-plane, with $+y$ being the North. “Top” and “Bottom” denote the $+z$ and $-z$ directions, respectively.

In addition to the regular electric and magnetic field samples, denoted by $E^{n}$ and $H^{n-(1/2)}$, respectively, we define the so-called hanging variables [15], $U^{n+(1/2)}$, which are magnetic field samples collocated with the electric fields on the boundary of the region. The hanging variables, which are tangential to the boundary and perpendicular to the corresponding electric field samples, are used to define the power supplied to the region and to facilitate the derivation of a self-contained dynamical model for the region based on FDTD update equations.

In this section, we take the FDTD equations for the regular field samples, $E^{n}$ and $H^{n-(1/2)}$, and cast them in the form of a dynamical system with the hanging variables, $U^{n+(1/2)}$, as its inputs. The output consists of electric field samples on the boundary of the region.

II-A Equations at Each Node

The temporal evolution of any electric field sample strictly inside the region, such as the $E_{x}|_{i+(1/2),j,k}^{n}$ sample shown in Fig. 2a, is described using a conventional FDTD update equation [2]

[TABLE]

where $\sigma_{x}$ and $\varepsilon_{x}$ are electrical conductivity and permittivity values, respectively, on the primary edge where the electric field is sampled. For clarity, we do not show the dependence of material properties on the location, although the proposed developments are valid in the most general case where all material properties are inhomogeneous. The iteration time step is denoted by $\Delta t$. Although $\Delta x$, the length of the edge where the sample is located, can be canceled in (1), we keep it for later derivations. The product $\Delta x\Delta y\Delta z$ is the volume associated with $E_{x}|_{i+(1/2),j,k}^{n}$, as illustrated in Fig. 3. We refer to the internal electric fields such as $E_{x}|_{i+(1/2),j,k}^{n}$ as electric fields of Type 1.

A conventional FDTD equation for the electric fields on the region’s faces, such as the South $E_{x}|_{i+(1/2),1,k}^{n}$ sample in Fig. 2b, would involve some magnetic field samples that are beyond the considered region. The use of the magnetic fields outside the region would require assumptions on the nature of the subsystems connected to the FDTD grid. Instead, we write a modified equation for the South electric field samples using the hanging variables on the South boundary

[TABLE]

Equation (2) is obtained from an FDTD-like approximation of Maxwell-Ampère law on a half-cell containing $E_{x}|_{i+(1/2),1,k}^{n}$. Using (2) instead of a conventional FDTD equation, we obtain a self-contained FDTD-like model for the fields in the region that does not involve field samples beyond its boundaries. This feature is crucial for investigating under which conditions the region is dissipative. It is also needed to obtain an FDTD model that can be connected to other subsystems, such as a grid of different resolution or a reduced model, where some magnetic field samples may not be available for a conventional FDTD equation. Samples like $E_{x}|_{i+(1/2),1,k}^{n}$ are classified as Type 2, which are the electric field samples on the faces of the boundary.

In a similar way, we write a modified FDTD equation for the samples located at the edges where two faces of the boundary come together. Those equations involve two hanging variables for each sample. For instance, in the case of the Bottom-South sample shown in Fig. 2c, the recurrence relation reads

[TABLE]

In (3), the $y$- and $z$-directed hanging variables, $U_{y}|_{i+(1/2),1,1}^{n+(1/2)}$ and $U_{z}|_{i+(1/2),1,1}^{n+(1/2)}$, respectively, are used to complete the line integral of magnetic field around $E_{x}|_{i+(1/2),1,1}^{n}$. The sample $E_{x}|_{i+(1/2),1,1}^{n}$ is an example of a Type 3 electric field, which is a field shared between two faces of the region’s boundary. The three types of $E_{y}^{n}$ and $E_{z}^{n}$ fields and their recurrence relations are defined analogously to $E_{x}^{n}$ and to (1)–(3).

The recurrence relation for the magnetic $H_{x}^{n-(1/2)}$ samples that are strictly inside the region is the conventional update FDTD equation, which reads

[TABLE]

where $\mu_{x}$ is the magnetic permeability at the sampling location. Equation for the boundary magnetic field sample $H_{x}|_{1,j+(1/2),k+(1/2)}^{n-(1/2)}$, which is located at a secondary edge of length $\Delta x/2$ normal to the West boundary, involves a factor of $\Delta x/2$, as opposed to $\Delta x$

[TABLE]

and similarly for the boundary faces on the East side. Recurrence relations for $H_{y}^{n-(1/2)}$ and $H_{z}^{n-(1/2)}$ are obtained in a similar fashion.

The coefficients $\Delta x\Delta y\Delta z$ in (4) and $(\Delta x/2)\Delta y\Delta z$ in (5) correspond to the dimensions of the volume associated with each $H_{x}^{n-(1/2)}$ sample, as shown in Fig. 3.

II-B Compact Matrix Form

In order to shorten the notation, we write the equations (1)–(5) for all field samples in a matrix form. Collecting all recurrence relations for $H_{x}^{n-(1/2)}$ samples (4) and (5) into matrix equations, we obtain

[TABLE]

where vectors $\mathbf{H}_{x}^{n-(1/2)}$ and $\mathbf{E}_{y}^{n}$ contain all $H_{x}^{n-(1/2)}$ and $E_{y}^{n}$ samples in the region, respectively. The coefficient matrix $\mathbf{D}_{\mu_{x}}$ is a diagonal matrix containing the magnetic permeability values on the edges where the corresponding elements of $\mathbf{H}_{x}^{n-(1/2)}$ are sampled. Matrix $\mathbf{G}_{zy}$, which consists of zeros, $+1$’s and $-1$’s, is a discrete $z$-derivative operator for $\mathbf{E}_{y}^{n}$. Similarly, $\mathbf{G}_{yz}$ is a discrete $y$-derivative operator for $\mathbf{E}_{z}^{n}$. Matrix $\mathbf{D}_{l^{\prime}_{x}}$ is a diagonal matrix containing the length of each secondary edge associated with samples in $\mathbf{H}_{x}^{n-(1/2)}$, namely $\Delta x$ for internal samples and $\Delta x/2$ for the samples on the East and West boundaries. Diagonal matrix $\mathbf{D}_{A_{x}}$ contains the area of primary faces where $H_{x}^{n-(1/2)}$ is sampled, which reduces to $\Delta y\Delta z\mathbf{I}_{N_{H_{x}}}$ for a uniform discretization, where $\mathbf{I}_{m}$ denotes an $m\times m$ identity matrix and $N_{H_{x}}$ is the number of $H_{x}^{n-(1/2)}$ samples in the region. As a result, the product $\mathbf{D}_{l^{\prime}_{x}}\mathbf{D}_{A_{x}}$ contains the volumes shown in the right panel of Fig. 3, associated with each sample in $\mathbf{H}_{x}^{n-(1/2)}$. Matrices $\mathbf{D}_{l_{y}}$ and $\mathbf{D}_{l_{z}}$ contain the length of the primary edges where the elements of $\mathbf{E}_{y}^{n}$ and $\mathbf{E}_{z}^{n}$ are sampled, respectively. These matrices reduce to $\Delta y\mathbf{I}_{N_{E_{y}}}$ and $\Delta z\mathbf{I}_{N_{E_{z}}}$ for uniform discretization.

In a similar fashion, we write in matrix form the recurrence relations for $\mathbf{H}_{y}^{n-(1/2)}$ and $\mathbf{H}_{z}^{n-(1/2)}$ as

[TABLE]

[TABLE]

with the coefficients defined similarly to those in (6).

From (1), (2), and (3), the recurrence relation for the $E_{x}^{n}$ samples can be written in matrix form as

[TABLE]

where $\mathbf{D}_{\varepsilon_{x}}$ and $\mathbf{D}_{\sigma_{x}}$ are diagonal matrices with the values of permittivity and electric conductivity, respectively, on the $x$-directed primary edges. The diagonal matrix $\mathbf{D}_{A^{\prime}_{x}}$ contains the area of the secondary faces where $E_{x}^{n}$ is sampled, which is equal to $\Delta y\Delta z$ for Type 1 nodes, $\Delta y\Delta z/2$ for Type 2 nodes, and $\Delta y\Delta z/4$ for Type 3 nodes. Matrix $\mathbf{B}_{{x,B}}$ consists of [math]’s and $1$’s that extract from the vector $\mathbf{U}_{y,B}^{n+(1/2)}$ the hanging variables corresponding to a Type 2 or Type 3 sample in $\mathbf{E}_{x}^{n}$. Matrices $\mathbf{B}_{{x,T}}$, $\mathbf{B}_{{x,S}}$, and $\mathbf{B}_{{x,N}}$ for the Top, South, and North faces serve an analogous purpose. In a similar manner, we write the matrix equations for $\mathbf{E}_{y}^{n}$ and $\mathbf{E}_{z}^{n}$.

Equations (6)–(8) describing the recurrence relation for all regular magnetic field variables can be written as

[TABLE]

where vectors $\mathbf{E}^{n}$ and $\mathbf{H}^{n-(1/2)}$ collect all electric and regular magnetic field variables, respectively

[TABLE]

The diagonal matrices $\mathbf{D}_{l^{\prime}}$ and $\mathbf{D}_{A}$ contain the secondary edge lengths and primary face areas associated with samples in $\mathbf{H}^{n-(1/2)}$. Matrix $\mathbf{D}_{\mu}$ contains permeability on the edges where samples in $\mathbf{H}^{n-(1/2)}$ are located. Coefficient matrix $\mathbf{C}$ is a discrete curl operator

[TABLE]

and $\mathbf{D}_{\l}$ is a diagonal matrix containing the length of the primary edges corresponding to the samples in $\mathbf{E}^{n}$.

Similarly, we obtain the following recurrence relation for the electric samples in the region by collecting (9) and similar equations for $\mathbf{E}_{y}^{n}$ and $\mathbf{E}_{z}^{n}$

[TABLE]

where $\mathbf{D}_{A^{\prime}}$ contains the secondary face areas associated with the samples in $\mathbf{E}^{n}$. Matrices $\mathbf{D}_{\varepsilon}$, and $\mathbf{D}_{\sigma}$ contain permittivity and electrical conductivity on the edges that correspond to the samples in $\mathbf{E}^{n}$. Matrix $\mathbf{Q}$ contains $1$’s and [math]’s that extract the hanging variables from vector $\mathbf{U}^{n+(1/2)}$. For Type 1 samples in $\mathbf{E}^{n}$, the corresponding rows in $\mathbf{Q}$ contain only zeros. For a Type 2 sample in $\mathbf{E}^{n}$, the row of $\mathbf{Q}$ contains a single $1$ in the column corresponding to the hanging variable that is involved in an equation such as (2) for that electric field sample. For a Type 3 sample, the row of $\mathbf{Q}$ contains two $1$’s in the columns corresponding to the two hanging variables in the equation like (3). Matrix $\mathbf{P}$ is a diagonal matrix that selects the signs for the hanging variables according to (9). The product $\mathbf{Q}\mathbf{P}$ has the following structure

[TABLE]

where

[TABLE]

Matrix $\mathbf{D}_{l^{\prime}_{U}}$ in (13) is a diagonal matrix containing the length of the edges where the hanging variables are sampled, namely $\Delta x$, $\Delta y$, or $\Delta z$ for the hanging variables collocated with Type 2 electric fields and $\Delta x/2$, $\Delta y/2$, or $\Delta z/2$ for the hanging variables associated with the electric fields of Type 3.

II-C Dynamical System Formulation

Matrix equations (10) and (13) can be written in the form of a dynamical system as follows

[TABLE]

where the state vector $\mathbf{x}^{n}$ contains the regular FDTD samples in the region

[TABLE]

The hanging variables $U^{n+(1/2)}$ are not included in the state vector, and instead are regarded as the input of the FDTD region

[TABLE]

where vector $\mathbf{U}_{E_{x}}^{n+(1/2)}$ collects the variables associated with the $E_{x}^{n}$ samples, namely the $U_{y}^{n+(1/2)}$ variables tangential to the Top and Bottom faces and $U_{z}^{n+(1/2)}$ variables for the North and South faces. Vectors $\mathbf{U}_{E_{y}}^{n+(1/2)}$ and $\mathbf{U}_{E_{z}}^{n+(1/2)}$ are defined similarly. The output vector $\mathbf{y}^{n}$ is defined as

[TABLE]

where $\mathbf{Y}_{E_{x}}^{n}$ contains the $E_{x}^{n}$ samples of Types 2 and 3 associated with the hanging variables in $\mathbf{U}_{E_{x}}^{n+(1/2)}$. Type 2 samples are included in $\mathbf{Y}_{E_{x}}^{n}$ once, while each Type 3 sample is included in $\mathbf{Y}_{E_{x}}^{n}$ twice, at the positions of the two corresponding hanging variables in $\mathbf{U}_{E_{x}}^{n+(1/2)}$. Vectors $\mathbf{Y}_{E_{y}}^{n}$ and $\mathbf{Y}_{E_{z}}^{n}$ are defined similarly to $\mathbf{Y}_{E_{x}}^{n}$.

Matrices $\mathbf{R}$, $\mathbf{F}$, $\mathbf{B}$, and $\mathbf{L}^{T}$ in (16a)–(16b) are given by

[TABLE]

The definitions given in this section allow us to write the FDTD equations for a 3-D region in the compact matrix form (16a)–(16b). Remarkably, equations (16a)–(16b) have the same structure as in the 2-D case [13]. This fact allows us to generalize previous results valid in two dimensions to the most general 3-D case.

III Dissipativity of a 3-D FDTD Region

In this section, we derive dissipativity conditions for the 3-D FDTD system (16a)–(16b) defined in Sec. II and explain their physical significance.

Definition 1**.**

According to the theory of dissipative systems [16], the discrete-time dynamical system (16a)–(16b) is dissipative with the supply rate $s\left(\mathbf{y}^{n},\mathbf{u}^{n+(1/2)}\right)$ if we can find a function $\mathcal{E}$ that satisfies

[TABLE]

The storage function $\mathcal{E}{\left(\mathbf{x}^{n}\right)}$ quantifies the energy stored in the FDTD region and the supply rate $s{\left(\mathbf{y}^{n},\mathbf{u}^{n+(1/2)}\right)}$ is the energy absorbed by the region through its boundaries between time points $n$ and $n\!+\!1$.

We define the storage function and the supply rate as

[TABLE]

[TABLE]

The diagonal matrix $\mathbf{S}$ in (23) is given by

[TABLE]

where $\mathbf{D}_{l_{Y}}$ contains the lengths of the primary edges on which the electric samples in $\mathbf{y}^{n}$ are taken. Expressions (22) and (23) generalize those from the 2-D case [13].

III-A Physical Meaning of the Supply Rate and Storage Function

By substituting (17) and (20a) into (22), we can write the storage function as

[TABLE]

which, using (10) to simplify the term inside the square brackets, reduces to

[TABLE]

This expression for the stored energy has been proposed in [11] for FDTD/FIT. Equation (26) reveals the physical meaning of the storage function (22) as a discrete analogy of

[TABLE]

which is the energy stored in electric and magnetic fields inside a given volume.

For the supply rate, we can substitute (24) into (23), obtaining

[TABLE]

The product $\mathbf{D}_{l_{Y}}\mathbf{D}_{l^{\prime}_{U}}$ is a diagonal matrix containing the area of boundary faces associated with each hanging variable and the corresponding electric field sample. Diagonal matrix P has a ${+1}$ on the faces where the Poynting vector points into the region and a ${-1}$ otherwise. As a result, (28) is a discrete version of the Poynting integral over the region’s boundary from time sample $n$ to $n\!+\!1$

[TABLE]

where $\hat{n}$ is the inward unit normal vector. Therefore, supply rate (23) is the total energy that enters the region from its boundaries between time $n$ and time $n\!+\!1$.

III-B Dissipativity Conditions for a 3-D FDTD Region

By substituting (22) and (23) into (21c), we arrive at the following theorem, which gives some simple conditions for the dissipativity of (16a)–(16b).

Theorem 1**.**

A system in the form (16a)–(16b) is dissipative with $\mathcal{E}{\left(\mathbf{x}^{n}\right)}$ in (22) as storage function and $s{\left(\mathbf{y}^{n},\mathbf{u}^{n+(1/2)}\right)}$ in (23) as supply rate if

[TABLE]

Proof.

See [13]. ∎

The physical meaning of the first condition (30a) can be better understood by applying the Schur complement [17] to transform (30a) into

[TABLE]

The first matrix inequality in (31) is true because lengths, areas, and permittivities in the region are positive. Applying a congruence with matrix $2(\Delta t\mathbf{D}_{l^{\prime}}\mathbf{D}_{A}\mathbf{D}_{\mu})^{-(1/2)}$ to the left hand side of the second inequality in (31), we obtain an equivalent condition to (31):

[TABLE]

where

[TABLE]

Condition (32) holds when

[TABLE]

where $s_{k}$ are the non-zero singular values of $\mathbf{\Sigma}$. Thus, condition (30a) can be seen as a generalized CFL limit, since it is applicable to the most general case of a lossy region with nonuniform material properties. Moreover, with a strategy similar to [11] we can show that a sufficient condition for (30a) to hold is that the classical FDTD CFL limit [2] is met

[TABLE]

where $\varepsilon$ is the smallest primary edge permittivity in the region and $\mu$ is the smallest secondary edge permeability. When the CFL limit is violated, the energy stored in some cells is no longer bounded below by zero, which can make them capable of supplying unlimited energy to the rest of the system, potentially causing instability.

Condition (30b) expands into

[TABLE]

which holds when the conductivity on each edge is non-negative.

The third condition (30c) is always true. In order to see this, we expand $\mathbf{L}\mathbf{S}$ as

[TABLE]

Right-multiplication by $\mathbf{D}_{l_{Y}}$ has the same effect on $\mathbf{Q}$ as left-multiplication by $\mathbf{D}_{\l}$. Hence, $\mathbf{Q}\mathbf{D}_{l_{Y}}=\mathbf{D}_{\l}\mathbf{Q}$, and as a result $\mathbf{L}\mathbf{S}=\mathbf{B}$.

In summary, Theorem (1) shows that the FDTD system (16a)–(16b) associated to an arbitrary region is dissipative if

all conductivities are non-negative, as one may obviously expect; 2. 2.

the CFL limit (30a) is respected.

If these two conditions are satisfied, the FDTD region can be arbitrarily interconnected with any other dissipative subsystem without violating stability.

IV Systematic Method for Stability Enforcement

The developments in the previous sections provide a powerful approach to construct both simple and advanced FDTD schemes with guaranteed stability. The method is general since it is applicable to complicated setups involving multiple subsystems, such as FDTD grids with different resolution, multiple boundary conditions, circuit models, or reduced-order models. The approach involves the following steps:

Hanging variables, which are seen as inputs, are defined at the boundary of each block. The corresponding electric field samples are regarded as outputs, as done in Sec. II for a 3D-FDTD region. 2. 2.

If subsystems cannot be connected directly, for example because of different grid resolution, an interpolation rule is created, and viewed as a separate subsystem. 3. 3.

Stability is enforced by ensuring that all subsystems are dissipative. For dynamical systems of the form (16a)–(16b) this can be done using (30a)–(30c). Since the connection of dissipative systems is also dissipative [18], the resulting scheme is stable by construction. 4. 4.

The most restrictive time step is taken to ensure that all subsystems are dissipative.

The proposed approach significantly simplifies the creation of new FDTD schemes with guaranteed stability, since dissipativity conditions can be imposed on individual subsystems. In contrast, existing techniques to analyze and enforce stability, such as the energy [11] and iteration [2, 12] methods, require the analysis of the entire scheme, which can be a formidable task. A remarkable feature of the proposed approach is its modularity. Given a set of subsystems that have been proven to be dissipative, any of their combination is guaranteed to be stable. This feature is hoped to accelerate scientific research in the FDTD area, since it allows researchers to create new stable schemes without having to redo stability analysis for every change. The proposed modular framework is relevant also for commercial FDTD solvers, where users want to be able to combine available FDTD subsystems in the way most suitable to model their problem, while having an absolute guarantee of stability.

It should be noted that, although we use a matrix formulation to investigate the dissipativity of FDTD equations, the final scheme can be implemented with scalar FDTD equations for optimal efficiency.

V Application to Stable 3-D FDTD Subgridding

As a demonstration of the stability framework in Sec. IV, we derive a stable subgridding algorithm for 3-D FDTD. In the proposed algorithm, a coarse grid is refined in selected regions to better resolve fine geometrical details. We denote the refinement factors as $r_{x}$, $r_{y}$, and $r_{z}$ in the $x$, $y$, and $z$ directions, respectively. The refinement factors are odd integers greater than one.

For simplicity, we present the proposed method by considering the scenario in Fig. 4, where the interface between the coarse and fine grids is planar. In the Appendix, we discuss how corners can be similarly treated.

In order to develop a stable subgridding scheme with the method in Sec. IV, we view a subgridding algorithm as a connection of four subsystems: boundary conditions, the coarse grid, the fine grid, and the interpolation rule that relates fields at the interface of the two grids, as depicted in Fig. 4.

V-A Interpolation Conditions

We first describe the interpolation rule established between the coarse and fine grid fields, which is the core of every subgridding algorithm. In the scenario of Fig. 4, the interpolation rule has to properly relate the fields tangential to the North boundary of the coarse grid to the fields tangential to the South boundary of the fine grid. We discuss the interpolation rule in detail for the $E_{z}^{n}$-$U_{x}^{n+(1/2)}$ pairs. The derivations for $E_{x}^{n}$-$U_{z}^{n+(1/2)}$ pairs can be done analogously.

As shown in the left panel of Fig. 5, we consider the portion of the interface surrounding one coarse electric field sample (denoted as $E_{z}^{n}$) and the corresponding hanging variable (denoted as $U_{x}^{n+(1/2)}$). The fine samples that fall in the same region in the case of $r_{x}\!=r_{z}\!=3$ are shown in the right panel of Fig. 5.

The fine samples are numbered form $1$ to $9$, as the nodes where they are sampled. The fine electric fields $\hat{E}_{z}^{n}$ are set equal in the $z$ direction

[TABLE]

and the fine hanging variables $\hat{U}_{x}^{n+(1/2)}$ are set equal in the $x$ direction

[TABLE]

where the “$\hat{\phantom{a}}$” notation refers to the variables of the refined region.

The coarse samples $E_{z}^{n}$ and $U_{x}^{n+(1/2)}$ are forced to be equal to the average of the fine samples through the following interpolation rules

[TABLE]

where

[TABLE]

and $\mathbf{T}_{m}$ is an $m\times 1$ matrix of ones.

Interpolation rules (38), (39), (40), and (41) draw inspiration from reciprocity theory for stable subgridding [19]. In the $z$ direction, the pair of rules (38) and (41) is similar to [19], except the magnetic fields are related directly at the interface, as opposed to a plane that is offset from the interface as in [19]. In the $x$ direction, the pair (39) and (40) is analogous, except the magnetic and electric interpolation rules are switched.

V-B Stability Proof

The stability of the proposed subgridding algorithm can be proved with the dissipation theory in Sec. IV. As shown in Fig. 4, the proposed scheme can be seen as the connection of various subsystems, which must all be dissipative in order to ensure stability. As proven in Sec. III, the coarse and fine grids are dissipative under their respective CFL limits. The perfect electric conductor (PEC) boundary condition is lossless, since imposing a null electric field means that no energy will be exchanged at any time.

What is left to investigate is whether the interpolation rule is dissipative or not. We show that the interpolation rule in Sec. V-A is actually lossless for any time step, by proving that the supply rate (28) of the interpolation rule is zero. The supply rate of the interpolation rule is the summation of terms, each one associated to one hanging variable on its boundaries. The terms associated with the $x$-directed hanging variables $U_{x}^{n+(1/2)}$ and $\hat{\mathbf{U}}_{x}^{n+(1/2)}$ read

[TABLE]

Symbol “$\otimes$” in (43) denotes the Kronecker product [20]. The supply rate term (43) accounts for the contribution of the $U_{x}^{n+(1/2)}$ and $\hat{U}_{x}^{n+(1/2)}$ samples in Fig. 5 and the discussion for $z$-directed hanging variables is analogous. Using the properties of the Kronecker product [20], we can rewrite (43) as

[TABLE]

Substituting (40) and (41) into (44), we see that the terms corresponding to the coarse and fine grids cancel each other

[TABLE]

Iterating the argument for all other coarse hanging variables and the corresponding groups of $r_{x}r_{z}$ hanging variables on the fine side, we conclude that the total supply rate to the interpolation rule is zero. Hence, the energy leaving the fine grid during a time step is equal to the energy that the interpolation rule supplies to the coarse grid. This makes the interpolation rule a lossless, and hence dissipative, subsystem. Since the coarse and fine grids are dissipative under their own CFL limits, the entire scheme is stable under the CFL limit of the fine grid, which is the most restrictive CFL limit.

V-C Practical Implementation

In this section we discuss how the proposed subgridding algorithm can be implemented. For updating all field samples strictly inside the two grids, one can use conventional FDTD equations. The update equation for the field samples on the refinement interface is instead derived from the interpolation rules (38), (39), (40), and (41).

As discussed in Sec. II-A, each electric field sample $E_{z}^{n}$ on the North boundary of the coarse grid satisfies

[TABLE]

where only those subscripts that are different from $i$, $j$, and $k\!+\!(1/2)$ are explicitly shown for clarity. Each electric field sample at a fine node ${\hat{m}\in\{1,2,\dots,9\}}$ with coordinates ${(\hat{\imath},\hat{\jmath}=\!1,\hat{k}\!+\!(1/2))}$ satisfies

[TABLE]

where a sample \hat{H}_{y\hat{m}}\big{|}^{n+(1/2)}_{\hat{\imath}+(1/2)} has the same coordinates as the node $\hat{m}$, except for the $x$ coordinate, which is $\hat{\imath}\!+\!(1/2)$. Similar notation is used for the other samples in (47).

Using (38a) and averaging (47) over the $\hat{k}$ indexes in the shaded area of Fig. 5, we obtain the following equation for the fine sample $\hat{E}_{z1}^{n}$

[TABLE]

and similarly for the other samples in $\hat{\mathbf{E}}_{z}^{n}$, namely for $\hat{E}^{n}_{z2}$ and $\hat{E}_{z3}^{n}$. Next, we use (41) to replace the average of the fine hanging variables in (48) with $U_{x}^{n+(1/2)}$

[TABLE]

Because of (39), the same substitution can be done in the equations for $\hat{E}^{n}_{z2}$ and $\hat{E}_{z3}^{n}$.

Writing (49) and the equations for $\hat{E}^{n}_{z2}$ and $\hat{E}_{z3}^{n}$ in matrix form, we obtain the following relation for the $\hat{\mathbf{E}}^{n}$ vector

[TABLE]

where

[TABLE]

and similarly for $\hat{\mathbf{H}}_{y,\hat{\imath}+(1/2)}^{n+(1/2)}$. Diagonal matrix $\hat{\mathbf{D}}_{\varepsilon_{z}}$ in (50) contains the permittivity of the half-cells adjacent to the fine grid’s South boundary, averaged in a way similar to the magnetic fields in (51a)–(51b). The conductivity matrix $\hat{\mathbf{D}}_{\sigma_{z}}$ is defined analogously.

Substituting (40) into (46) and multiplying the result on the left by $\mathbf{T}_{r_{x}}/r_{x}$, we obtain

[TABLE]

where $\mathbf{1}_{{r_{x}}}$ denotes a square matrix of ones of size ${{r_{x}}\!\times{r_{x}}}$. Adding (50) to (52) to eliminate $U_{x}^{n+(1/2)}$, we obtain an ${{r_{x}}\!\times{r_{x}}}$ matrix equation

[TABLE]

with the vector of unknowns $\hat{\mathbf{E}}_{z}^{n+1}$, where matrices $\mathbf{M}_{\varepsilon_{z}}$ and $\mathbf{M}_{\sigma_{z}}$ are given by

[TABLE]

Equation (53) relates the electric field vector $\hat{\mathbf{E}}_{z}$ to the surrounding magnetic fields in both the coarse and the fine grids. We can see that this last step eliminates the hanging variables on the boundaries of the two grids, which serve as a temporary means to connect the two grids. Equation (53) is a recursive relation for the electric field samples at the interface, that can be used to compute $\hat{\mathbf{E}}_{z}^{n+1}$ knowing the electric and magnetic samples at the previous times. Although (53) is not fully explicit in terms of $\hat{\mathbf{E}}_{z}^{n+1}$, the matrix in front of this unknown is very small, having size ${r_{x}\times r_{x}}$. Moreover, it is a constant. Therefore, this matrix has to be inverted only once, leading to an explicit update equation for $\hat{\mathbf{E}}_{z}^{n+1}$, which will be used for all time steps.

Overall, the proposed subgridding scheme consists of the following steps.

Starting from $E^{n}$ and $H^{n+(1/2)}$, we update all electric fields strictly inside fine and coarse regions using standard FDTD equations such as (1). 2. 2.

We update the electric fields on the interface between the two grids using (53) and a similar equation for $\mathbf{E}_{x}^{n}$ samples. Once (53) is solved, the fine electric samples $\hat{E}_{z4}^{n+1}$, $\ldots$ , $\hat{E}_{z9}^{n+1}$ can be found using the equalities (38a)–(38c) and the coarse electric field sample can be updated using (40). 3. 3.

All regular magnetic field samples are updated using standard FDTD equations, such as (4) and (5), to obtain the values at $n\!+\!(3/2)$.

The inclusion of material properties throughout the derivations rigorously guarantees stability in the case when arbitrary permittivities and conductivities are assigned to cells on each side of the subgridding interface. In contrast, the stability of many existing schemes with material traverse can only be verified numerically [21, 22] and some schemes exhibit inaccuracy when objects intersect the interface [23].

VI Numerical Examples

The subgridding algorithm from Sec. V was implemented in MATLAB in order to test the validity of the proposed theory for FDTD stability and to assess the accuracy of the proposed subgridding method. Time-consuming portions of the code, such as the update equations for the fields inside each grid, were written using vectorized operations.

VI-A Stability Verification

We verify the stability of the proposed algorithm by simulating the subgridding scenario shown in Fig. 6 for a million time steps. The setup consists of a 12 cm $\!\times\!$ 12 cm \!\times\!\leavevmode\nobreak\12 cm cavity with PEC walls. The cavity is discretized with a uniform coarse mesh with $\Delta x\!=\!\Delta y\!=\!\Delta z\!=\!$ 1 cm. A 4 cm $\!\times\!$ 4 cm $\!\times\!$ 4 cm subgridding region with refinement ratio $r\!=\!5$ in all directions is placed at the center of the cavity. In order to test the correctness of stability analysis in the case when material properties vary, we randomly assign permittivity values to cells in both grids between the free space value $\varepsilon_{0}$ and $3\varepsilon_{0}$ using the $rand()$ function in MATLAB. In another simulation, we test if losses are properly handled by introducing, in addition to the varying permittivities, random conductivities between zero and 5$\times$10*-5* S/m. We excite the cavity using a Gaussian pulse with half-width at half-maximum bandwidth of 3.53 GHz. Time step is set to 99% of the fine grid’s CFL limit.

From the results, presented in Fig. 7, it can be seen that no signs of instability occur during the million time steps of the simulation, even if nonuniform materials traverse the coarse-fine interface. This result validates the correctness of the stability analysis framework presented in this paper.

VI-B Material Traverse

In this simulation, we test whether the proposed subgridding method produces consistent results when the subgridding interface traverses material boundaries, which is known to be a challenging scenario for subgridding schemes [23]. We use the setup shown in Fig. 8, where the subgridding region is placed in three different ways: in Case 1 and Case 3, a block of material traverses different sides of the interface. In the reference Case 2, the block is fully enclosed by the subgridding boundary. We perform the test for a copper block and for a lossy dielectric block ($\varepsilon_{r}\!=\!3$, $\sigma\!=\!0.05$ S/m). A 15-cm perfectly matched layer (PML) terminates the domain. The coarse grid is set to $\Delta x\!=\!\Delta y\!=\!\Delta z\!=\!$ 1 cm and the subgridding region is refined by a factor of $r\!=\!3$ in each direction. A Gaussian source with half-width at half-maximum bandwidth of 1.02 GHz is used. The iteration time step is 99% of the CFL limit of the fine grid.

Fig. 9 shows the electric field recorded at the probe in the three cases. The two cases with material traverse are in very good agreement with the reference simulation, showing that the proposed algorithm does not lose accuracy when different materials traverse the subgridding boundary. The test is successful for both lossy dielectric and for copper.

VI-C Meta-Screen

In this example, we use subgridding to study the three-slot meta-screen shown in Fig. 10. The setup for this simulation is similar to [24], except for differences in discretization and screen material. In particular, we set the plate conductivity to 1.3$\times$106 S/m and thickness of the plate to 6 mils, based on the value in [25]. The slot in the center is 13.2 mm by 1.2 mm and the two satellite slots on the sides are 17 mm by 0.6 mm, placed at a center-to-center distance of 3 mm away from the central slot. These dimensions were designed to achieve subwavelength focusing at 10 GHz [25]. We place a line of probes 4.572 mm away from the meta-screen, which corresponds to approximately 15% of the wavelength at 10 GHz – distance at which measurements were performed in [25]. The simulation region was 25.2 mm$\times$19.2024 mm$\times$43.4 mm including PMLs.

As in [24], we terminate the region before the meta-screen with PEC walls on the East and the West sides, and with perfect magnetic conductor walls on the Top and the Bottom sides. A five-cell PML layer is added on all sides in order to mimic an unbounded domain.

The satellite slots have widths of only 2% of the wavelength, which makes the problem intrinsically multiscale. Moreover, [24] reports resonance effects that make it difficult to resolve the structure near the design frequency. We choose this example to investigate the performance of the proposed subgridding scheme.

We set the coarse grid to $\Delta x\!=\!\Delta z\!=\!$ 0.7 mm, and $\Delta y\!=\!$ 0.1524 mm. In order to properly resolve the meta-screen structure and the fields, we refine a 9.8 mm$\times$2.5908 mm$\times$21 mm area around the slots with $r_{x}\!=\!r_{z}\!=\!7$ and $r_{y}\!=3$, using the proposed subgridding method. As a reference, we perform a conventional FDTD simulation with a nonuniform grid, where a region of dimensions 9.8 mm$\times$1.0668 mm$\times$21 mm is refined by $r_{x}\!=\!r_{z}\!=\!7$ and $r_{y}\!=3$ by means of varying $\Delta x$, $\Delta y$, and $\Delta z$. We also simulate the case where the entire structure is discretized with the coarse mesh.

A modulated Gaussian waveform with 10 GHz central frequency and 8.24 GHz half-width at half-maximum bandwidth is used to excite the structure with a sheet of $x$-directed uniform current density on the South side. The proposed method and FDTD with nonuniform discretization are run at $\Delta t=$ 0.1362 ps, whereas the coarse grid simulation has the iteration time step of 0.4810 ps. These time steps correspond to 99% of the CFL limit of each mesh. The waveform is recorded for 5.5 ns, as in [24].

Fig. 11 shows the amplitude of the $x$ component of electric field at 10 GHz, which was obtained with a Fourier transform of the time-domain waveforms. The results from subgridding are in very good agreement with the reference simulation. In contrast, the coarse resolution was insufficient to perform an accurate simulation of the meta-screen. This is not surprising, since any mesh with $\Delta x$ larger than 0.3 mm or $\Delta z$ larger than 0.1 mm cannot correctly resolve the meta-screen structure, let alone the fields around the slots.

Simulation times, as well as the number of primary FDTD cells are shown in Table I. Subgridding almost doubles the speed-up compared to the nonuniform grid refinement, which is already much faster than the simulation with uniform grid. This is a result of the lower number of unknowns associated with subgridding. The guarantee of stability of the proposed method allows for the reduction of simulation time without sacrificing reliability. Moreover, the ability of the method to handle material traverse allows placing the subgridding region across the meta-screen.

VII Conclusion

We proposed a dissipation theory for 3-D FDTD inspired by the theory of dissipative dynamical systems [14]. The theory shows that the FDTD equations for the fields in a 3-D region can be interpreted as a dynamical system that is dissipative when time step satisfies the CFL limit.

The proposed dissipation theory provides a powerful and systematic method to create new FDTD schemes with guaranteed stability. By showing that each component of the scheme is dissipative (such as grids, boundary conditions, or embedded lumped components), the overall coupled scheme is guaranteed to be dissipative, and thus stable. This approach makes stability analysis simpler, more intuitive, and modular, since dissipation conditions are imposed on each component independently and are based on the familiar concept of energy.

With the new framework, we develop a stable 3-D subgridding scheme supporting material traverse. Numerical examples show good performance of the scheme when accelerating multiscale problems, as well as its ability to produce accurate and stable simulations when various materials traverse the subgridding interface. The dissipation framework is not restricted to subgridding algorithms, and could be applied to developing more advanced schemes in the future. As a proof of concept, we have applied the dissipation framework in 2-D to guarantee the stability of FDTD coupled to reduced-order models [18].

In this section, we handle the case when the fine grid is surrounded by the coarse grid on all six sides. In this scenario, the locations where two faces of the boundary come together require a special treatment. Fig. 12 shows the portion of the fine grid’s boundary on its South-West edge, which is discussed in detail in this section.

-A Interpolation Conditions

Unlike in Sec. V, the edge case involves two sets of hanging variables: those in the $y$ direction on the West side and those in the $x$ direction on the South side. As in Sec. V, we interpolate the fields on the shaded portion of the interface so that the $z$-directed electric field samples are equal in the $z$ direction

[TABLE]

The fine $y$-directed hanging variables are set equal in the $y$ direction and the fine $x$-directed hanging variables are set equal in the $x$ direction

[TABLE]

We collect the distinct electric field samples into vectors as follows

[TABLE]

where vector $\hat{\mathbf{E}}_{z(W)}^{n}$ contains the distinct samples on the West half of the shaded region in Fig. 12, except for the South-West sample at ${(1,1,\hat{k}\!+\!(1/2))}$ and similarly for $\hat{\mathbf{E}}_{z(S)}^{n}$ on the South side. Variable $\hat{E}_{z(SW)}^{n}$ is always the single sample at ${(1,1,\hat{k}\!+\!(1/2))}$, which is excluded from $\hat{\mathbf{E}}_{z(W)}^{n}$ and $\hat{\mathbf{E}}_{z(S)}^{n}$. In general, $\hat{\mathbf{E}}_{z(W)}^{n}$ and $\hat{\mathbf{E}}_{z(S)}^{n}$ are vectors of size ${(r_{y}\!-\!1)/2}$ and ${(r_{x}\!-\!1)/2}$, respectively, although in the $r_{x}\!=\!r_{y}\!=3$ example they reduce to scalar quantities.

Similarly, the distinct hanging variables in (57a)–(57c) are collected into vectors $\hat{\mathbf{U}}_{y}^{n+(1/2)}$ and $\hat{\mathbf{U}}_{x}^{n+(1/2)}$

[TABLE]

[TABLE]

Analogously to (40), we force the coarse electric field sample $E_{z}^{n}$ to be equal to the average of $\hat{E}_{z}^{n}$ over the West portion of the shaded area in Fig. 12

[TABLE]

with the weighting coefficients chosen to correctly account for the area of boundary faces associated with each $\hat{E}_{z}^{n}$ sample. The column vector of ones $\mathbf{T}_{W}$ has size ${(r_{y}-1)/2}$. A similar condition is applied on the South portion of the shaded area in Fig. 12

[TABLE]

As in the planar case, the coarse hanging variables are forced to be equal to the average of the corresponding fine hanging variables

[TABLE]

-B Stability Proof

In this section, we show that with the interpolation rule in Sec. -A for the edges of the boundary, the interpolation rule remains a lossless subsystem. Consider the portion of the coarse-fine interface consisting of the South shaded area in Fig. 12 on the fine side and of the adjacent $\Delta x/2\times\Delta z$ boundary face on the coarse side. For this portion of the interface, the net contribution of the $x$-directed hanging variables to the supply rate from the fine and coarse grids to the interpolation rule boils down to

[TABLE]

analogously to (44). Substituting (64), (65) becomes

[TABLE]

which is zero because of (62). The proof for other hanging variables near the edges of the boundary can be done with a similar argument. This shows that even in presence of corners, the interpolation rule is a lossless system, since its supply rate is zero.

-C Practical Implementation

The update equations for the electric fields in the shaded area in Fig. 12 are derived using the modified FDTD equations on each side of the interface and the interpolation conditions (56a)–(56c), (57a)–(57c), (61), (62), (63), and (64). The coarse $E_{z}^{n}$ field at ${(i,j,k\!+\!(1/2))}$ needs to satisfy the modified FDTD equation

[TABLE]

Applying a similar averaging procedure as in Sec. V-C, we obtain the following relations for the fine electric field vectors defined in (58)

[TABLE]

where

[TABLE]

where $(\hat{\imath}\!=\!1,\hat{\jmath}\!=\!1,\hat{k}\!+\!(1/2))$ is the coordinate of node $\hat{m}$. In general, vectors $\hat{\mathbf{H}}_{x(W),\hat{\jmath}\pm(1/2)}^{n+\frac{1}{2}}$ and $\hat{\mathbf{H}}_{y(W),\hat{\imath}+(1/2)}^{n+(1/2)}$ have size ${(r_{y}\!-\!1)/2}$, with each element corresponding to a sample in $\hat{\mathbf{E}}_{z(W)}$. Quantities $\hat{H}_{y(SW),\hat{\imath}+(1/2)}^{n+(1/2)}$, $\hat{\mathbf{H}}_{y(S),\hat{\imath}\pm(1/2)}^{n+(1/2)}$, and $\hat{\mathbf{H}}_{y(W),\hat{\imath}+(1/2)}^{n+(1/2)}$ are defined similarly to (69a)–(69d). Permittivity and conductivity coefficients $\hat{\varepsilon}_{z(SW)}$, $\hat{\sigma}_{z(SW)}$, $\hat{\mathbf{D}}_{\varepsilon_{z(W)}}$, $\hat{\mathbf{D}}_{\sigma_{z(W)}}$, $\hat{\mathbf{D}}_{\varepsilon_{z(S)}}$, and $\hat{\mathbf{D}}_{\sigma_{z(S)}}$ are obtained from averaging the values in the half-cells adjacent to the boundary over $\hat{k}$ indexes. Equations (68a)–(68c) serve a purpose similar to (50).

We now combine (67) for the $z$-directed electric field on the coarse side, the recurrence relations (68a)–(68c) for the fields on the fine side, and the interpolation rules (61) and (62), into a linear system of the form

[TABLE]

with an ${(r_{x}/2+r_{y}/2+3)\times 1}$ vector of unknowns

[TABLE]

Solving (70) for $\mathbf{z}^{n+1}$ gives values of $E_{z}^{n+1}$, $\hat{E}_{z(SW)}^{n+1}$, $\hat{\mathbf{E}}_{z(W)}^{n+1}$, and $\hat{\mathbf{E}}_{z(S)}^{n+1}$. With (56a)–(56c), samples $\hat{E}_{z3}^{n+1}$ and $\hat{E}_{z5}^{n+1}$ are obtained from $\hat{E}_{z1}^{n+1}$, samples $\hat{E}_{z4}^{n+1}$ and $\hat{E}_{z6}^{n+1}$ from $\hat{E}_{z2}^{n+1}$, and $\hat{E}_{z4^{\prime}}^{n+1}$ and $\hat{E}_{z6^{\prime}}^{n+1}$ from $\hat{E}_{z2^{\prime}}^{n+1}$. The size of system (70) is only $(r_{x}/2+r_{y}/2+3)\times(r_{x}/2+r_{y}/2+3)$ and the resulting overhead is typically small.