1. Introduction

The inverse design method, which is a combined simulation and optimization technique, has achieved remarkable results as regards the design of high-performance nanophotonic devices [1–6]. In this process, the processing speed and result accuracy are directly dependent on the simulation method. For nanophotonic device design using this approach, the finite-difference time-domain (FDTD) method or the finite difference frequency domain method (FDFD) can be employed to simulate the electromagnetic field [1,7]. Compared with the FDTD method, the FDFD method can be extremely fast and efficient, with lower memory consumption. In addition, the FDFD method exhibits superior numerical dispersion analysis and better ability to import oblique incidences and resolve sharp resonances [8–10]. As a result of these many advantages, the FDFD method is primarily employed for electromagnetic simulations in the context of inverse design.

In inverse design, optimization can be performed using the adjoint variable method, through which the gradient of an objective function for any arbitrary number of design parameters can be calculated via two electromagnetic simulations [11]. In this approach, a level-set function is usually used to design a more practical, discretely valued nanophotonic device [12–14]. When the level-set function is incorporated, adjoint optimization is usually referred to as topology optimization. Adjoint/topology optimization has been applied to linear, nonlinear, and active devices, and has exhibited significant ability as regards the design of nanophotonic devices at optical wavelengths [15–17].

In the inverse design process, the design time increases significantly with the design area; therefore, design of large-scale devices within a limited timeframe is difficult. Various studies have been conducted to reduce the design process time consumption for an enlarged design area, thereby accelerating the process. For example, researchers from Stanford University previously optimized the parameters of a perfectly matched layer (PML) [9] to improve the convergence of the iterative solver for an FDFD simulation. They then applied a 2.5-dimension FDFD method to the inverse design process to accelerate this process for a photonic crystal fiber [8]. They also employed the Schur complement method [18,19] to reduce the coefficient matrix size, and the conjugate orthogonal conjugate gradient method to reduce the iteration time of each step [20].

However, inverse design of large-scale devices remains challenging, as the design area is important and must be simulated at high spatial resolution. The background area containing waveguides and the PML can be simulated at low spatial resolution. As different spatial resolutions can be applied to different areas, an adaptive mesh is more efficient than a structured mesh.

Non-uniform gridding, octree grids, and subgridding/adaptive meshes, hereafter collectively referred to as “adaptive meshes,” are usually used to accelerate electromagnetic simulation. Several studies considering the FDFD method and adaptive meshes have been conducted to date [21–24]. For example, Kulas and Mrozowski proposed a novel three-dimensional subgridding scheme applicable to the finite-difference technique in the time and frequency domains [21]. When the refinement factor was set to 3, the field values on the coarse mesh overlapped with some of the field values on the fine mesh; hence, mapping of the field values from the coarse to fine mesh was no longer necessary. Therefore, the reflection between coarse and fine meshes could be reduced. Subsequently, Podwalski, Kowalczyk, and Mrozowski used a subgridding scheme to improve the efficiency of FDFD simulations for microwave and radiofrequency (RF) devices [22]. Further, Fotyga et al. presented a novel local and nested (multilevel) macromodel that could further improve the FDFD solver speed [23], while Luan et al. proposed a novel FDFD method for RF device simulation based on octree grids, combining a linear interpolation scheme and a second-order interpolation scheme to reduce reflection [24]. In those works, the simulation was accelerated through introduction of an adaptive mesh. However, the interpolation scheme on the boundary between the fine and coarse meshes introduced numerical dispersion and artificial reflections to the simulation results. In a small number of studies, the FDFD method has also been applied to an adaptive mesh for optical device simulation [23]. To the best of the authors’ knowledge, however, incorporation of an adaptive mesh with the FDFD method has not yet been applied to inverse design.

To apply FDFD simulation to an adaptive mesh, an interpolation scheme between the coarse and fine meshes must be implemented on their boundary. Linear, high-order, and moving least-squares meshless interpolation have been used in previous studies [24–27]. However, the moving least-squares meshless interpolation, in which weights are applied to different simulation types, is not suitable for inverse design. Further, the simulation results obtained using linear interpolation are not sufficiently accurate. High-order interpolation yields greater accuracy; however, the simulation speed is also reduced.

In this study, to accelerate the inverse design process for large-scale devices, an adaptive mesh is introduced to an FDFD system and optimization is realized based on the adjoint variable method. As regards the interpolation scheme, experiments conducted in this study show that the second-order interpolation scheme is sufficiently accurate and, hence, a second-order interpolation scheme is employed. To implement the inverse design process, a linear coefficient matrix based on the FDFD method must be constructed. A less efficient approach to obtaining the coefficient matrix involves addition of several specific points to the coefficient matrix on the boundary between the coarse and fine meshes. However, the present study provides a general framework that allows effective construction of a coefficient matrix on an adaptive mesh. For a 1 × 2 beam splitter at 1550 nm as an inverse design target, a numerical experiment is performed and the FDFD simulation accuracy is verified through comparison with the FDTD module of the Lumerical commercial software package. Furthermore, to assess the speed of the design process, an FDFD simulation with a structured mesh is conducted for comparison.

The remainder of this paper is organized as follows. Section 2.1 presents the basic theory of the adjoint variable method and Section 2.2 illustrates the interpolation scheme. In Section 2.3, the general framework for constructing the FDFD coefficient matrix is presented. Section 3 analyzes the results of the inverse design numerical experiment and the limitations of the inverse design process based on the adaptive mesh. Finally, Section 4 concludes this study.

2. Basic theory and principles

2.1 FDFD method and adjoint variable method

Before discussing acceleration of the inverse design process, we first briefly describe the FDFD and adjoint variable methods, which are used to realize the inverse design process.

The FDFD method based on the frequency-domain Maxwell equations is the simulation method typically applied for performance analysis of nanophotonic devices in inverse design. For this method, simulation is performed by solving a linear equation with form

Equation (1) can be simplified as

In the inverse design process, the permittivity distribution with respect to the electric field can also be determined. This process can usually be realized using a convex optimization algorithm, with the adjoint variable method being used to calculate the optimization gradient.

For the latter, the figure of merit, described by a scalar objective function L, must be maximized for the design variables φ. The adjoint variable method then computes the sensitivity. Here,

Typically, ∂L/∂E can be solved analytically. In Eq. (3), φ represents the material distribution in the design region. To determine ∂E/∂φ, Eq. (2) can be derived as follows:

Hence, the following can be obtained:

Substituting Eq. (5) into Eq. (3) yields

In each optimization step, the design variables must be updated from φ_old to φ_new as follows:

2.2 Interpolation scheme

An adaptive mesh is usually a combination of several meshes having different spatial resolutions. In this study, only two meshes are used to construct the adaptive mesh, having low and high spatial resolutions and referred to as the coarse and fine meshes, respectively. Distributions of structured and adaptive meshes are illustrated in Fig. 1. Hereafter, the Yee cell sizes of the coarse and fine meshes are labeled dc and df, respectively.

 

Fig. 1. Mesh distributions of a) structured and b) adaptive meshes.

Download Full Size | PDF

In inverse design, the designed device is usually a thin slice with fewer Yee cells in the z-direction. Thus, it is not necessary to set different spatial resolutions for the coarse and fine meshes in the z-direction, and only the interpolation scheme for the other two directions is discussed in this section. To expand the application potential of the proposed approach, the interpolation scheme in three directions is presented in Appendix 1. To perform FDFD simulation on an adaptive mesh, field values are obtained from the fine mesh with the field values from the coarse mesh, and field values are obtained from the coarse mesh with the field values from the fine mesh.

First, we discuss acquisition of field values from the fine mesh to the coarse mesh. In the FDFD simulation, the field component is mapped onto a staggered mesh, which is usually referred to as the Yee mesh. To apply an adaptive mesh to an FDFD system, the coarse and fine meshes must obey the Yee theory, which states that different field components in a single Yee cell have different positions. Therefore, if one entire Yee cell is split into multiple sub-Yee cells, an interpolation scheme is required to directly obtain the field values from the sub-Yee cells for the single entire Yee cell. In the context of the present work, if the entire Yee cell is regarded as the Yee cell on the coarse mesh and the multiple sub-Yee cells are taken as the Yee cells on the fine mesh, the field values must be interpolated from the fine mesh to the coarse mesh. However, previous research has shown that each field value on the coarse mesh can be spatially aligned with one of the field values on the fine mesh, as shown in Fig. 2, if the refinement factor (dc/df) is set to 3 [21]. Then, the field values on the coarse mesh can be obtained by determining the spatially aligned field values on the fine mesh. Based on this theory, the interpolation scheme can be simplified and the reflection on the boundary between meshes can be reduced.

 

Fig. 2. Electric (E) and magnetic (H) field arrangement on adaptive mesh for refinement factor of 3.

Download Full Size | PDF

Next, interpolation to the fine mesh with field values from the coarse mesh is considered. To simplify the explanation, only the interpolation scheme for the Hz field component is discussed, where H represents the magnetic field. Using a similar approach, interpolation functions can be constructed for the Ex, Ey, Ez, Hx, and Hy components. The interpolation function for Hz in the Yee cell can be defined as

 

Fig. 3. Spatial arrangement of five Hz components in two dimensions.

Download Full Size | PDF

To calculate a, b, c, d, and e, it is assumed that Hz can be expressed as a second-order function of the spatial coordinates such that

By substituting the field values and spatial positions of the five Hz components into Eq. (10), the following can be obtained:

Then, solving for a’, b’, c’, d’, and e’ and substituting the solutions into Eq. (10), the interpolation function can be written as

Compared with Eq. (9), a, b, c, d, and e can be expressed as

Based on the interpolation function for each component, several sub-interpolation matrices can be constructed for each component, here labeled D_Ex, D_Ey, D_Ez, D_Hx, D_Hy, and D_Hz. Details of the sub-interpolation matrix construction are provided in Appendix 2. By combining D_Ex, D_Ey, D_Ez, D_Hx, D_Hy, and D_Hz, interpolation matrices labeled D_inp, D_E, and D_H can be obtained, such that

It must be emphasized that, in the present study, all the field values are mapped from the coarse mesh to the fine mesh to simplify the implementation. However, only parts of D_inp are used on the boundary to construct the FDFD system, such that the interpolation scheme is only applied on the boundary.

2.3 Implementation

In this section, we discuss construction of an FDFD system based on the interpolation scheme. To simplify the construction, two FDFD systems are built: one on the coarse mesh and the other on the fine mesh. We define Ce_coarse and Ch_coarse, and Ce_fine and Ch_fine as the difference matrices of the two FDFD systems, respectively. To simplify the discussion, the region of the design area in which the Yee cells on the fine mesh spatially align with the Yee cells on the coarse mesh is defined as the overlapping area. The regions of the design area containing the other Yee cells that are not spatially aligned are defined as the non-overlapping areas. The spatial distribution of the overlapping areas, non-overlapping areas, and background area of the fine mesh is shown in Fig. 4.

 

Fig. 4. Spatial distributions of Yee cells on a) coarse and b) fine meshes.

Download Full Size | PDF

Based on this spatial distribution, Ce_fine, Ce_coarse, and D_E can be expressed as

To build the FDFD system on an adaptive mesh, the following matrices should be constructed: two difference matrices representing the curl operator on the electric and magnetic fields, Ce_combine and Ch_combine, respectively; two diagonal matrices representing the permittivity and permeability, ε_combine and μ_combine, respectively; and one vector representing the electric current source density, J_combine. Of all the required matrices and vectors, construction of Ce_combine and Ch_combine is particularly challenging.

First, we consider construction of the difference matrix. Only construction of Ce_combine is discussed in detail, as Ch_combine can be constructed using a similar method. To simplify the construction, Ce_combine is obtained by combining Ce_coarse, Ce_fine, and D_E. This combination can be performed because the design area has higher spatial resolution than the background area in the inverse design process. Therefore, only the portions of the difference matrix related to the background area on the coarse mesh and the design area on the fine mesh are considered.

For detailed discussion of the difference matrix construction, we briefly review the interpolation scheme. The Hz distribution shown in Fig. 5 is used as an example. On the coarse mesh, the three field components marked with hollow circles are required in order to realize the forward difference process. As shown in Fig. 5(a), the hollow-circle field components overlap some of the field components on the fine mesh. In this case, the hollow-circle field components can be obtained by selecting the overlapping components on the fine mesh. On the fine mesh, the nine field components marked with hollow circles are required in order to realize the backward difference process. As shown in Fig. 5(b), the hollow-circle field components can be obtained by implementing an interpolation scheme.

 

Fig. 5. Interpolation schemes from a) fine to coarse mesh and b) from coarse to fine mesh.

Download Full Size | PDF

First, a scheme for interpolation from the fine to coarse mesh is considered. Because the Yee cells in the overlapping area on the fine mesh are spatially aligned with those in the design area on the coarse mesh, the coupling matrix Cec_d→b can be used directly. Combining Cec_background, Cec_d→b, and Ce_fine yields the following:

Next, interpolation from a coarse to fine mesh is considered. To introduce the interpolation scheme, D_E is incorporated into Ce_add1 to yield

Further, ε_combine, μ_combine, and J_combine can be constructed by extracting and combining components of the permittivity, permeability, and electric current source from the FDFD systems on the coarse and fine meshes having the spatial distribution shown in Fig. 4.

Using the matrices Ce_combine, Ch_combine, ε_combine, and μ_combine, the coefficient matrix for the FDFD system on the adaptive mesh can be constructed as

The field value of the adaptive mesh can be obtained by solving a linear system in the form

Following introduction of the interpolation matrix to the FDFD system, the coefficient matrix is no longer the Hermitian matrix. However, for the experiment conducted in this study, the biconjugate gradient (BICG), which is usually applied to Hermitian system solution, could provide an appropriate solution for the FDFD system on the adaptive mesh. Therefore, the BICG solver was chosen for application to the FDFD system.

Finally, as shown in Eq. (20), ε_combine, which incorporates the design variable φ, is on the diagonal line of A_combine. Therefore, A_combine can be directly introduced to the adjoint variable method.

3. Numerical experiment

This section describes the numerical experiment and discusses the results. All computations were performed on a computer with two graphics processing units (NVIDIA Tesla P100 16G) and one central processing unit (Intel Xeon E5-2630). The software was based on the g++ platform with the CUDA package.

Because inverse design acceleration is more advantageous for large-scale problems, only an FDFD system containing more than 180 × 180 × 50 Yee cells was considered. Thus, the Yee cell dimensions were set to 75 and 25 nm for the coarse and fine meshes, respectively, where the dimensions were identical in both the x- and y-directions; the z-direction dimension for both meshes was set to 25 nm. On the structured mesh, the Yee cell dimensions in both the x- and y-directions were set to 25 nm. Therefore, all simulation areas for the inverse design process exceeded 4.5 μm × 4.5 μm × 1.25 μm.

To illustrate the advantage of the adaptive mesh, a 1 × 2 beam splitter was used as the experimental target. The initial scheme of the beam splitter is illustrated in Fig. 6. The waveguide width and Si layer thickness were 500 and 220 nm, respectively. In the inverse design process, the Si and SiO₂ refractive indices were set to 3.5 and 1.5, respectively. To normalize the comparison, the distance between the design area and the boundary of the simulation area was set to 1.125 μm.

 

Fig. 6. Initial scheme of 1 × 2 beam splitter. a) Top view and b) Front view.

Download Full Size | PDF

3.1 Acceleration

In this study, the simulation area dimensions in both the x- and y-directions were used as the length fraction. Before analyzing the acceleration, the size of the coefficient matrix for the adaptive-mesh FDFD system and was compared with that for the structured-mesh FDFD system under different length fractions, as shown in Fig. 7. The coefficient matrix size for the adaptive-mesh FDFD system was far smaller than that of the structured-mesh FDFD system. As the length fraction increased, the difference in coefficient matrix size between the structured- and adaptive-mesh FDFD systems increased. For the BICG solver, the speed of each step depends on the coefficient matrix size; thus, the above result implies that, for each step of the BICG solver, the calculation for the adaptive-mesh FDFD system can be faster than that for the structured mesh.

 

Fig. 7. Coefficient matrix size comparison between adaptive- and structured-mesh FDFD systems.

Download Full Size | PDF

For the BICG solver, the convergence also directly affects the computational time. Here, the number of iterative steps required for convergence was taken to represent convergence. A comparison of the convergence for different length fractions is shown in Fig. 8, and reveals that the iterative solver convergence on the adaptive mesh was poorer than that on the structured mesh for a simulation area smaller than 6.1 μm × 6.1 μm × 1.25 μm. As the simulation area increased, the iterative-solver convergence degraded for both cases. However, for a simulation area exceeding 6.1 μm × 6.1 μm × 1.25 μm, the iterative-solver convergence on the adaptive mesh exceeded that on the structured mesh.

 

Fig. 8. Convergence comparison between adaptive- and structured-mesh FDFD systems.

Download Full Size | PDF

Each iterative-solver step was conducted more quickly on the adaptive mesh than on the structured mesh; thus, the BICG solver and, hence, the inverse design process, can be implemented more quickly on an adaptive mesh for large-scale problems.

To analyze the performance of the adaptive-mesh FDFD simulation for a large-scale problem, the total computational times and matrix-vector-multiplication times of both the adaptive- and structured-mesh FDFD simulations were compared for different length fractions, as shown in Fig. 9. Even though the simulation area was quite large, the adaptive-mesh FDFD simulation was faster. Furthermore, the comparative results for the total computational time and vector multiplication time had a similar form, which confirms that the iterative-solver time consumption depends on the vector multiplication time. As the vector multiplication time depends on the coefficient matrix size, the time-consumption comparison shown in Fig. 9 matches the size comparison shown in Fig. 7.

 

Fig. 9. Time-consumption comparison between adaptive- and structured-mesh FDFD simulations. a) Matrix vector multiplication time and b) total BICG-solver computational time for one simulation process or gradient calculation process.

Download Full Size | PDF

3.2 Accuracy verification

To verify the accuracy of the simulation results, the FDTD module of the Lumerical commercial optical simulation software was used to simulate a nanophotonic device, which was designed with an inverse design process on an adaptive mesh. As shown in Fig. 10, the simulation result obtained from the adaptive-mesh FDFD method matched that of the Lumerical FDTD module. However, for several steps, the transmission obtained with the adaptive mesh differed slightly from the Lumerical result. This phenomenon arose because the interpolation scheme provided an approximate result only, and because reflection was introduced on the boundary between the coarse and fine meshes. The design and simulation results are presented in Fig. 11.

 

Fig. 10. Transmission comparison between adaptive-mesh FDFD simulation and Lumerical FDTD simulation for inverse design process. a) Top and b) bottom ports.

Download Full Size | PDF

 

Fig. 11. Simulation result for nanophotonic devices designed with Lumerical FDTD module on structured mesh.

Download Full Size | PDF

Usually, improved FDFD simulation accuracy corresponds to reduced speed, and vice versa. Here, the accuracy and acceleration were compared to appropriately analyze the FDFD simulation efficiency. Table 1 compares the transmission in the top and bottom ports for the designed 1 × 2 beam splitters and the total time consumption of the inverse design process with different mesh resolutions. Note that the transmission was calculated using the Lumerical FDTD module. For condition normalization, the target coupling efficiency at both output ports was set to 0.49–0.51. The maximum optimization step was defined as 30. The inverse design process with the adaptive mesh achieved a 1.79-fold acceleration over the structured-mesh case (25 nm). In addition, for a mesh resolution of approximately 40 nm, the time consumption (i.e., speed) of the inverse design process was similar to that of the adaptive mesh. However, the performance of the nanophotonic device obtained through the inverse design process on the structured mesh was lower than that obtained through the inverse design process on the adaptive mesh.

Table 1. Comparison of design-result transmissions for different mesh resolutions

View Table | View all tables in this article

3.3 Discussion

Inverse design on an adaptive mesh can accelerate the design process and increase the design area. However, the adaptive mesh introduces a small reflection at the boundary between the coarse and fine meshes, which generates a difference in transmission. To verify the simulation results, transmissions in the top and bottom ports given by the adaptive-mesh FDFD simulation and Lumerical FDTD simulation were compared for different frequencies. As shown in Fig. 12, both the resonant frequency and transmission differed slightly between the two simulations. This phenomenon arose because the interpolation scheme introduced a reflection at the boundary between the coarse and fine meshes. Further, the simulations with different spatial resolutions may have had different resonant frequencies. This problem can be solved by employing a high-order interpolation scheme.

 

Fig. 12. Beam-splitter transmission comparison between adaptive-mesh FDFD simulation and Lumerical FDTD simulation for different frequencies. a) Top and b) bottom ports.

Download Full Size | PDF

For inverse design, the designed device performance degrades as the spatial resolution decreases. To analyze the limitations of the adaptive mesh with regard to the spatial resolution, the designed-device transmissions for different spatial resolutions were compared. As detailed in Table 2, when the coarse mesh resolution was 105 nm, the top- and bottom-port transmissions were 0.345 and 0.347, respectively; these values are low for the design of a nanophotonic device. However, for a 90-nm coarse mesh resolution, the top- and bottom-port transmissions were acceptable, at 0.413 and 0.414, respectively. Therefore, we suggest that the coarse mesh resolution must exceed 90 nm for effective implementation of the inverse design process.

Table 2. Transmission comparison for different spatial resolutions

View Table | View all tables in this article

4. Conclusion

In this study, the inverse design process computation was accelerated using an FDFD simulation on an adaptive mesh. Experiment results showed that 1.79-fold acceleration was achieved by the adaptive-mesh FDFD simulation over the structured-mesh case for an FDFD system simulation area of 4.5 μm × 4.5 μm × 1.25 μm. For verification, the simulation results were compared with those given by the Lumerical FDTD module, with good agreement being obtained. Further, a 2.25 μm × 2.25 μm beam splitter at 1550 nm was obtained through inverse design using an FDFD simulation on an adaptive mesh, with 0.458-dB insertion loss and low time consumption.

However, FDFD simulation on an adaptive mesh usually yields a system error at the boundary between the fine and coarse meshes, which in turn causes artificial reflection and numerical dispersion in the inverse design process. Therefore, in the context of inverse design, we believe that an adaptive mesh is more suitable for initial high-precision designs of nanophotonic devices rather than final designs.

In addition to acceleration and simulation area expansion, FDFD simulation on an adaptive mesh can be used to develop machine-learning-based optimization methods to identify the more important parts of a nanophotonic device. Hence, components requiring individual simulations with high and low resolutions could be determined. The acceleration achieved in this work could then be further improved by optimizing the device on multilevel meshes. In addition, inverse design on an adaptive mesh could be used to determine the most important characteristic in the design area, for a high-resolution, regularly shaped area. Indeed, a regular shape could be adopted as a characteristic for nanophotonic devices. The authors intend to develop these techniques in future work.

Appendix 1. Interpolation scheme in three directions

The interpolation scheme in three directions differs from that in two directions. First, instead of 3 × 3 field values on the fine mesh, there are 3 × 3 × 3 field values. Second, the interpolation function has three inputs. The interpolation function can be written as

(22)$$\begin{aligned} &Hz({x,y,z} )= a \times Hz({{x_0},{y_0} - 1,{z_0}} )+ b \times Hz({{x_0} - 1,{y_0},{z_0}} )+ c \times Hz({{x_0},{y_0},{z_0}} )+ \\ &d \times Hz({{x_0},{y_0},{z_0} - 1} )+ e \times Hz({{x_0},{y_0},{z_0} + 1} )+ f \times Hz({{x_0} + 1,{y_0},{z_0}} )+ g \times Hz({{x_0},{y_0} + 1,{z_0}} )\end{aligned}. $$

The spatial distribution of Hz(x, y, z) is shown in Fig. 13. To calculate a, b, c, d, e, f, and g, spatial coordinates are used to represent Hz, where

(23)$$Hz({x,y,z} )= a^{\prime}\Delta {x^2} + b^{\prime}\Delta {y^2} + c^{\prime}\Delta {z^2} + d^{\prime}\Delta x + e^{\prime}\Delta y + f^{\prime}\Delta z + g^{\prime}, $$

and

(24)$$\begin{aligned}\Delta x &= x - {x_0}\\ \Delta y &= y - {y_0}\\ \Delta z &= z - {z_0} \end{aligned}. $$

The interpolation function must satisfy

(25)$$\left\{ {\begin{array}{l} {Hz({{x_0},{y_0} - 1,{z_0}} )= b^{\prime} - e^{\prime} + g^{\prime}}\\ {Hz({{x_0} - 1,{y_0},{z_0}} )= a^{\prime} - d^{\prime} + g^{\prime}}\\ {Hz({{x_0},{y_0},{z_0}} )= g^{\prime}}\\ {Hz({{x_0},{y_0},{z_0} - 1} )= c^{\prime} - f^{\prime} + g^{\prime}}\\ {Hz({{x_0},{y_0},{z_0} + 1} )= c^{\prime} + f^{\prime} + g^{\prime}}\\ {Hz({{x_0} + 1,{y_0},{z_0}} )= a^{\prime} + d^{\prime} + g^{\prime}}\\ {Hz({{x_0},{y_0} + 1,{z_0}} )= b^{\prime} + e^{\prime} + g^{\prime}} \end{array}} \right.. $$

 

Fig. 13. Interpolation scheme for Hz component in three directions.

Download Full Size | PDF

From Eq. (25), a, b, c, d, e, f, and g can be calculated such that

(26)$$\left\{ {\begin{array}{l} {a = \left( {\frac{{\Delta {y^2}}}{2} - \frac{{\Delta y}}{2}} \right)}\\ {b = \left( {\frac{{\Delta {x^2}}}{2} - \frac{{\Delta x}}{2}} \right)}\\ {c = ({1 - \Delta {x^2} - \Delta {y^2} - \Delta {z^2}} )}\\ {d = \left( {\frac{{\Delta {z^2}}}{2} - \frac{{\Delta z}}{2}} \right)}\\ {e = \left( {\frac{{\Delta {z^2}}}{2} + \frac{{\Delta z}}{2}} \right)}\\ {f = \left( {\frac{{\Delta {x^2}}}{2} + \frac{{\Delta x}}{2}} \right)}\\ {g = \left( {\frac{{\Delta {y^2}}}{2} + \frac{{\Delta y}}{2}} \right)} \end{array}} \right.$$

Appendix 2. Construction of sub-interpolation matrix

Based on Eqs. (9) and (14), the interpolation matrix for a Yee cell on a coarse mesh can be constructed as

(27)$$DH{z_{one\_cell}} = \left[ {\begin{array}{ccccc} {\left( {\frac{{{y_1}^2}}{2} - \frac{{{y_1}}}{2}} \right)}&{\left( {\frac{{{x_1}^2}}{2} - \frac{{{x_1}}}{2}} \right)}&{({1 - {x_1}^2 - {y_1}^2} )}&{\left( {\frac{{{x_1}^2}}{2} + \frac{{{x_1}}}{2}} \right)}&{\left( {\frac{{{y_1}^2}}{2} + \frac{{{y_1}}}{2}} \right)}\\ \vdots & \vdots & \vdots & \vdots & \vdots \\ {\left( {\frac{{{y_9}^2}}{2} - \frac{{{y_9}}}{2}} \right)}&{\left( {\frac{{{x_9}^2}}{2} - \frac{{{x_9}}}{2}} \right)}&{({1 - {x_9}^2 - {y_9}^2} )}&{\left( {\frac{{{x_9}^2}}{2} + \frac{{{x_9}}}{2}} \right)}&{\left( {\frac{{{y_9}^2}}{2} + \frac{{{y_9}}}{2}} \right)} \end{array}} \right], $$

where x₁, …x₉ and y₁, …y₉ are the spatial coordinates of each Yee cell on the fine mesh within the Yee cell on the coarse mesh.

The interpolation matrix DH_{zone_cell} is the same for each Yee cell on a coarse mesh. Therefore, the sub-interpolation matrix DHz can be constructed through point-by-point setting of DH_{zone_cell}.

To simplify the discussion, the simulation area shown in Fig. 14 is taken as an example. Then, DHz can be expressed as

(28)

 

Fig. 14. Index of each Yee cell on coarse and fine meshes.

Download Full Size | PDF

To accelerate the construction, the indexes of each Yee cell on the fine and coarse meshes can be constructed. Faster construction of DHz can be achieved based on these indexes and DHz_{one_cell}. Further, the sub-interpolation matrices DEx, DEy, DEz, DHx, and DHy can be constructed based on a similar method.

Funding

National Natural Science Foundation of China (61935003).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz, T. M. Babinec, and J. Vučković, “Inverse design and demonstration of a compact and broadband on-chip wavelength demultiplexer,” Nat. Photonics 9(6), 374–377 (2015). [CrossRef]  

2. Z. Ye, J. Qiu, C. Meng, L. Zheng, Z. Dong, and J. Wu, “Inverse design of a SOI T-junction polarization Beamsplitter,” J. Phys.: Conf. Ser. 844, 012009 (2017). [CrossRef]  

3. N. V. Sapra, D. Vercruysse, L. Su, K. Y. Yang, J. Skarda, A. Y. Piggott, and J. Vuckovic, “Inverse design and demonstration of broadband grating couplers,” IEEE J. Select. Topics Quantum Electron. 25(3), 1–7 (2019). [CrossRef]  

4. N. Mohammadi Estakhri, B. Edwards, and N. Engheta, “Inverse-designed metastructures that solve equations,” Science 363(6433), 1333–1338 (2019). [CrossRef]  

5. B. Kailkhura, B. Gallagher, S. Kim, A. Hiszpanski, and T. Y. Han, “Reliable and explainable machine-learning methods for accelerated material discovery,” npj Comput. Mater. 5(1), 108–109 (2019). [CrossRef]  

6. L. Su, R. Trivedi, N. V. Sapra, A. Y. Piggott, D. Vercruysse, and J. Vučković, “Fully-automated optimization of grating couplers,” Opt. Express 26(4), 4023–4034 (2018). [CrossRef]  

7. B. Shen, P. Wang, R. Polson, and R. Menon, “An integrated-nanophotonics polarization beamsplitter with 2.4 × 2.4 μm2 footprint,” Nat. Photonics 9(6), 378–382 (2015). [CrossRef]  

8. J. Lu, S. Boyd, and J. Vučković, “Inverse design of a three-dimensional nanophotonic resonator,” Opt. Express 19(11), 10563–10570 (2011). [CrossRef]  

9. W. Shin and S. Fan, “Choice of the perfectly matched layer boundary condition for frequency-domain Maxwell’s equations solvers,” J. Comput. Phys. 231(8), 3406–3431 (2012). [CrossRef]  

10. W. Shin and S. Fan, “Accelerated solution of the frequency-domain Maxwell’s equations by engineering the eigenvalue distribution of the operator,” Opt. Express 21(19), 22578–22595 (2013). [CrossRef]  

11. M. G. Damigos and E. M. Papoutsis, “kgIachagias, Giannakoglou K C, adjoint variable-based shape optimization with bounding surface constraints,” Int. J. Numer. Meth. Fluids (2020).

12. N. Lebbe, ““Contribution in topological optimization and application to nanophotonics,” Analysis of PDEs,” NNT 2019GREAM047, tel-02518286 (2019).

13. R. E. Christiansen, F. Wang, and O. Sigmund, “Topological insulators by topology optimization,” Phys. Rev. Lett. 122(23), 234502 (2019). [CrossRef]  

14. S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vucković, and A. W. Rodriguez, “Inverse design in nanophotonics,” Nat. Photonics 12(11), 659–670 (2018). [CrossRef]  

15. T. W. Hughes, M. Minkov, I. A. D. Williamson, and S. Fan, “Adjoint method and inverse design for nonlinear nanophotonic devices,” ACS Photonics 5(12), 4781–4787 (2018). [CrossRef]  

16. J. Wang, Y. Shi, T. Hughes, Z. Zhao, and S. Fan, “Adjoint-based optimization of active nanophotonic devices,” Opt. Express 26(3), 3236–3248 (2018). [CrossRef]  

17. L. Su, A. Y. Piggott, N. V. Sapra, J. Petykiewicz, and J. Vučković, “Inverse design and demonstration of a compact on-chip narrowband three-channel wavelength demultiplexer,” ACS Photonics 5(2), 301–305 (2018). [CrossRef]  

18. N. Z. Zhao, S. Boutami, and S. Fan, “Accelerating adjoint variable method based photonic optimization with Schur complement domain decomposition,” Opt. Express 27(15), 20711–20719 (2019). [CrossRef]  

19. N. Zhao, S. Verweij, W. Shin, and S. Fan, “Accelerating convergence of an iterative solution of finite difference frequency domain problems via Schur complement domain decomposition,” Opt. Express 26(13), 16925–16939 (2018). [CrossRef]  

20. L. Su, D. Vercruysse, J. Skarda, N. V. Sapra, J. A. Petykiewicz, and J. Vučković, “Nanophotonic inverse design with SPINS: Software architecture and practical considerations,” Appl. Phys. Rev. 7(1), 011407 (2020). [CrossRef]  

21. L. Kulas and M. Mrozowski, “Low-reflection subgridding,” IEEE Trans. Microwave Theory Tech. 53(5), 1587–1592 (2005). [CrossRef]  

22. J. Podwalski, P. Kowalczyk, and M. Mrozowski, “Efficient multiscale finite difference frequency domain analysis using multiple macromodels with compressed boundaries,” Prog. Electromagn. Res. 126(463), 479 (2012). [CrossRef]  

23. G. Fotyga, P. Kowalczyk, L. Kulas, K. Nyka, J. Podwalski, and M. Mrozowski, “Reduced order models in computational electromagnetics (in memory of Ruediger Vahldieck),” ASIA PAC. INTERNATIONAL SYMPOSIUM ON ELECTROMAGNETIC COMPATIBILITY (APEMC), 2012.

24. H. Luan, Y. Geng, Y. Yu, and S. Guan, “Three-dimensional transient electromagnetic numerical simulation using FDFD based on octree grids,” IEEE Access 7, 161052–161063 (2019). [CrossRef]  

25. X. Zhu, Y. Shu, C. Luo, F. Huo, and W. Zhu, “Transient electromagnetic voltage imaging of dense UXO-like targets based on improved mathematical morphology,” IEEE Access 8, 150341–150349 (2020). [CrossRef]  

26. R. Goldade, Y. Wang, M. Aanjaneya, and C. Batty, “An adaptive variational finite difference framework for efficient symmetric octree viscosity,” ACM Trans. Graph. 38(4), 1–14 (2019). [CrossRef]  

27. A. Fedeli, C. Montecucco, and G. L. Gragnani, “Open-source software for electromagnetic scattering simulation: The case of antenna design,” Electronics 8(12), 1506 (2019). [CrossRef]