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

by   Fadime Bekmambetova, et al.

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.



page 1

page 2

page 3

page 4


On some stable boundary closures of finite difference schemes for the transport equation

We explore in this article the possibilities and limitations of the so-c...

Discrete transparent boundary conditions for the two-dimensional leap-frog scheme

We develop a general strategy in order to implement (approximate) discre...

Boundary treatment of implicit-explicit Runge-Kutta method for hyperbolic systems with source terms

In this paper, we develop a high order finite difference boundary treatm...

Stability analysis of a singlerate and multirate predictor-corrector scheme for overlapping grids

We use matrix stability analysis for a singlerate and multirate predicto...

Analysis of the SBP-SAT Stabilization for Finite Element Methods Part II: Entropy Stability

In the research community, there exists the strong belief that a continu...

A SBP-SAT FDTD Subgridding Method Using Staggered Yee's Grids Without Modifying Field Components

A summation-by-parts simultaneous approximation term (SBP-SAT) finite-di...

A Stable FDTD Subgridding Scheme with SBP-SAT for Transient Electromagnetic Analysis

We proposed a provably stable FDTD subgridding method for accurate and e...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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

Fig. 1: Illustration of the 3-D FDTD region with the hanging variables shown in green and the regular electric and magnetic samples in blue and red, respectively. The secondary grid is shown with the dashed line.

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,  [2]. 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 [15], , 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

(a) a(b) a(c)
a ab ac
Fig. 2: Variables involved in the recurrence relations for (a) internal samples (Type 1), (b) samples of Type 2 on the South face, and (c) the samples of Type 3 on the Bottom-South edge of the boundary. The common subscript “” is omitted for clarity.
Fig. 3: Illustration of the volumes (shown in dashed lines), associated with each electric (left panel) and magnetic (right panel) field sample in (1)–(5).

The temporal evolution of any electric field sample strictly inside the region, such as the sample shown in Fig. 2a, is described using a conventional FDTD update equation [2]


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.

The coefficients in (4) and in (5) correspond to the dimensions of the volume associated with each 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 samples (4) and (5) into matrix equations, we obtain


where vectors

and 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).

From (1), (2), and (3), the recurrence relation for the samples can be written in matrix form as


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 .

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


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

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


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 .

Matrices , , , and in (16a)–(16b) are given by


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 if we can find a function that satisfies


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


where contains the lengths of the primary edges on which the electric samples in 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


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


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


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


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

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 in (22) as storage function and in (23) as supply rate if


See [13]. ∎

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


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 [11] we can show that a sufficient condition for (30a) to hold is that the classical FDTD CFL limit [2] 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 .

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

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

  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:

  1. 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. 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. 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. 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 , , 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.

Fig. 4: Illustration of the subgridding scenario. Left: cross section of the primary grid. Right: subsystem interpretation.

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