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  is stable when the iteration time step is below the Courant-Friedrichs-Lewy (CFL) limit . 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  or circuit simulators . 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 , the iteration method  and the energy method . 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, 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 . The method is based on the theory of dissipative dynamical systems . 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 , , and primary cells in the , , and directions, respectively. Electric field is sampled on the primary grid edges at the integer time points, , and magnetic field is sampled on the primary grid faces at the time points shifted by half of a time step, . For simplicity, we assume that the region is discretized uniformly with primary cells of dimensions . As shown in Fig. 1, we use cardinal directions in the -plane, with being the North. “Top” and “Bottom” denote the and directions, respectively.
In addition to the regular electric and magnetic field samples, denoted by and , respectively, we define the so-called hanging variables , , 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, and , and cast them in the form of a dynamical system with the hanging variables, , as its inputs. The output consists of electric field samples on the boundary of the region.
Ii-a Equations at Each Node
where and 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 . Although , the length of the edge where the sample is located, can be canceled in (1), we keep it for later derivations. The product is the volume associated with , as illustrated in Fig. 3. We refer to the internal electric fields such as as electric fields of Type 1.
A conventional FDTD equation for the electric fields on the region’s faces, such as the South 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
Equation (2) is obtained from an FDTD-like approximation of Maxwell-Ampère law on a half-cell containing . 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
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
In (3), the - and -directed hanging variables, and , respectively, are used to complete the line integral of magnetic field around . The sample 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 and fields and their recurrence relations are defined analogously to and to (1)–(3).
The recurrence relation for the magnetic samples that are strictly inside the region is the conventional update FDTD equation, which reads
where is the magnetic permeability at the sampling location. Equation for the boundary magnetic field sample , which is located at a secondary edge of length normal to the West boundary, involves a factor of , as opposed to
and similarly for the boundary faces on the East side. Recurrence relations for and are obtained in a similar fashion.
Ii-B Compact Matrix Form
where vectorsand contain all and samples in the region, respectively. The coefficient matrix is a diagonal matrix containing the magnetic permeability values on the edges where the corresponding elements of are sampled. Matrix , which consists of zeros, ’s and ’s, is a discrete -derivative operator for . Similarly, is a discrete -derivative operator for . Matrix is a diagonal matrix containing the length of each secondary edge associated with samples in , namely for internal samples and for the samples on the East and West boundaries. Diagonal matrix contains the area of primary faces where is sampled, which reduces to for a uniform discretization, where denotes an identity matrix and is the number of samples in the region. As a result, the product contains the volumes shown in the right panel of Fig. 3, associated with each sample in . Matrices and contain the length of the primary edges where the elements of and are sampled, respectively. These matrices reduce to and for uniform discretization.
In a similar fashion, we write in matrix form the recurrence relations for and as
with the coefficients defined similarly to those in (6).
where and are diagonal matrices with the values of permittivity and electric conductivity, respectively, on the -directed primary edges. The diagonal matrix contains the area of the secondary faces where is sampled, which is equal to for Type 1 nodes, for Type 2 nodes, and for Type 3 nodes. Matrix consists of ’s and ’s that extract from the vector the hanging variables corresponding to a Type 2 or Type 3 sample in . Matrices , , and for the Top, South, and North faces serve an analogous purpose. In a similar manner, we write the matrix equations for and .
where vectors and collect all electric and regular magnetic field variables, respectively
The diagonal matrices and contain the secondary edge lengths and primary face areas associated with samples in . Matrix contains permeability on the edges where samples in are located. Coefficient matrix is a discrete curl operator
and is a diagonal matrix containing the length of the primary edges corresponding to the samples in .
Similarly, we obtain the following recurrence relation for the electric samples in the region by collecting (9) and similar equations for and
where contains the secondary face areas associated with the samples in . Matrices , and contain permittivity and electrical conductivity on the edges that correspond to the samples in . Matrix contains ’s and ’s that extract the hanging variables from vector . For Type 1 samples in , the corresponding rows in contain only zeros. For a Type 2 sample in , the row of contains a single 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 contains two ’s in the columns corresponding to the two hanging variables in the equation like (3). Matrix is a diagonal matrix that selects the signs for the hanging variables according to (9). The product has the following structure
Matrix in (13) is a diagonal matrix containing the length of the edges where the hanging variables are sampled, namely , , or for the hanging variables collocated with Type 2 electric fields and , , or for the hanging variables associated with the electric fields of Type 3.
Ii-C Dynamical System Formulation
where the state vector contains the regular FDTD samples in the region
The hanging variables are not included in the state vector, and instead are regarded as the input of the FDTD region
where vector collects the variables associated with the samples, namely the variables tangential to the Top and Bottom faces and variables for the North and South faces. Vectors and are defined similarly. The output vector is defined as
where contains the samples of Types 2 and 3 associated with the hanging variables in . Type 2 samples are included in once, while each Type 3 sample is included in twice, at the positions of the two corresponding hanging variables in . Vectors and are defined similarly to .
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 . 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
The storage function quantifies the energy stored in the FDTD region and the supply rate is the energy absorbed by the region through its boundaries between time points and .
We define the storage function and the supply rate as
The diagonal matrix in (23) is given by
Iii-a Physical Meaning of the Supply Rate and Storage Function
which, using (10) to simplify the term inside the square brackets, reduces to
which is the energy stored in electric and magnetic fields inside a given volume.
The product 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 on the faces where the Poynting vector points into the region and a otherwise. As a result, (28) is a discrete version of the Poynting integral over the region’s boundary from time sample to
where is the inward unit normal vector. Therefore, supply rate (23) is the total energy that enters the region from its boundaries between time and time .
Iii-B Dissipativity Conditions for a 3-D FDTD Region
See . ∎
The first matrix inequality in (31) is true because lengths, areas, and permittivities in the region are positive. Applying a congruence with matrix to the left hand side of the second inequality in (31), we obtain an equivalent condition to (31):
Condition (32) holds when
are the non-zero singular values of. 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  we can show that a sufficient condition for (30a) to hold is that the classical FDTD CFL limit  is met
where is the smallest primary edge permittivity in the region and 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
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 as
Right-multiplication by has the same effect on as left-multiplication by . Hence, , and as a result .
all conductivities are non-negative, as one may obviously expect;
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.
If subsystems cannot be connected directly, for example because of different grid resolution, an interpolation rule is created, and viewed as a separate subsystem.
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  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 , , and in the , , and 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 - pairs. The derivations for - 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 ) and the corresponding hanging variable (denoted as ). The fine samples that fall in the same region in the case of are shown in the right panel of Fig. 5.
The fine samples are numbered form to , as the nodes where they are sampled. The fine electric fields are set equal in the direction
and the fine hanging variables are set equal in the direction
where the “” notation refers to the variables of the refined region.
The coarse samples