Open Access
Add to BibSonomyAbstract
In this paper a new technique for numerical analysis of microstructured optical fibers is proposed. The technique uses a combination of model order reduction method and discrete function expansion technique. A significant reduction of the problem size is achieved (by about 85%), which results in much faster simulations (up to 16 times) without affecting the accuracy.
©2011 Optical Society of America
1. Introduction
This work is a continuation and extension of the approach presented in our previous papers [1, 2]. The new idea proposed herein is the application of model order reduction (MOR) technique to guided wave propagation problems. Model order reduction has originally been developed for acceleration of numerical analysis of dynamical systems for which one evaluates the transfer function by inverting a system matrix (or solving a system of equations) at many frequency points [3]. The general idea is to reduce the complexity of the problem (i.e. decrease the number of states which translates in a smaller matrix to be inverted at each frequency point) while preserving (within a limited frequency band) its input-output behavior. Among various techniques that have been developed for linear time invariant systems one can mention projection based methods such as the PVL (Pade via Lanczos) [4], PRIMA (Passive Reduced-Order Interconnect Macromodeling Algorithm) [5] and ENOR (Efficient Nodal Order Reduction) [6]. They all generate a small set of orthogonal vectors which span the solution space in the limited frequency range. Model order reduction can also be used in computational electromagnetics or photonics [7, 8, 9, 10, 11, 12]. It is particulary straightforward to apply this technique in the methods which use spatial discretization to solve Maxwell’s equations. These equations can be cast in the state-space form with state variables being the field samples at the discretization points. MOR allows one to eliminate internal redundant variables, and thus it is particularly useful when the discretization is very fine.
The MOR technique for dicretized Maxwell’s equations on Yee’s mesh can be implemented in three main steps [9]. The schematic idea of this process is presented in Fig. 1. In the first step a particular region of the domain is selected and its input/output ports are defined. In the FD method the whole domain is discretized and evenly covered with Yee’s mesh and the ports mentioned above can be associated with the field samples at the boundary of the region. Mathematically, the global matrix operator, formed from Maxwell’s grid equations, is separated into one part corresponding to the model and the second part corresponding to the rest of the domain. The second stage is a projection of the operator (its part corresponding to the model) onto a fixed “smaller” subspace. The projection eliminates most of internal state variables, while the relation between the output and input ports remains preserved for a limited frequency range. The projection basis is obtained automatically with ENOR algorithm [6] and it can by applied to the analysis of inhomogeneous regions. Finally, the reduced part of the operator is coupled with the unprojected part of the domain. As a result, the problem size is smaller and the computation time is much shorter. MOR has been successfully applied to reduce the complexity of FDTD (finite difference in time domain) and FDFD (finite difference in frequency domain) methods [9].
In this paper we present the application of MOR to analysis of propagation problems. It is important to note, that this type of problems leads to a linear space-invariant system and the space-state representation involves a set of differential equations involving space rather than the time variable. Consequently, the propagation constant, rather than frequency plays a role of a ”sweep” parameter. For this reason a different formulation of the system of equations to be reduced is required (different from the standard one considered in [13, 14, 1]).
An efficiency of the reduction strongly depends on a number of the ports coupling the model with the global domain (field samples at the boundary of the model). Therefore to improve the reduction efficiency a combination of MOR and discrete function expansion (DFE) technique is proposed. In this approach samples of the field at the boundary are replaced with a smaller number of harmonic amplitudes of the field. This improvement results in a significantly smaller size of the problem and in a very fast analysis (up to 16 times shorter than in the standard FD method).
2. Formulation
All structures analyzed in this paper are homogeneous along the z direction, so we can assume that E(ρ, ϕ,z) = E(ρ, ϕ)e −γz and H(ρ, ϕ,z) = H(ρ, ϕ)e −γz, where γ is the propagation coefficient. In the analysis a standard Yee’s mesh in the cylindrical coordinate system is applied (Fig. 2). We assume that the computational domain is discretised into M and N points in the radial and angular direction, respectively. Consequently, Yee’s mesh measures K = M × N points.
Following these assumptions, Maxwell’s curl equations can be expressed in a discrete matrix form, as [1]
where and Δρ, Δϕ are the space discretization in the ρ and ϕ direction, respectively.The elimination of longitudinal components of the fields from the equations (1) gives
where To invert the matrices in (9) no special procedure is required ( , are diagonal and , are orthogonal so their inverse is a transposed matrix).The formulation represented by equations (8) differs significantly from the standard one presented in [13, 14, 1]. However, also in this approach the eigenvalue problem for the propagation coefficient can be obtained (from a simple substitution of relations (8)). The main advantage of this formulation is that taking the inverse Laplace transform with respect to γ one gets a system of state-space equations which have the form that makes the application of the MOR technique possible.
In the first step a region for which MOR is to be applied must be selected (modeled area). At the boundary of that region a set of input and output ports must be defined. Fig. 3 shows the scheme of the coupling. The modeled area is denoted by a yellow background, the regular grid by white background and the coupling ports are enveloped with a thick line.
Fig. 3 Boundary between the model (yellow background) and a regular grid (white background) with coupling ports.
It is convenient to segregate the fields samples located inside the model (e M, h M) from the ones located outside (e U, h U). In such a case equations (8) can be rewritten in the following form
where , are the discrete operator for the grid point inside the model, , for outside and S ( e ), S ( h ) are coupling matrices [9].In order to apply a MOR technique one has to consider a relation between the fields inside and outside the model, which can be obtained from a combination of the equations (10)
The relation can be written in the form of a matrix-valued transfer function involving the boundary fields only where h B represents the field samples placed at the boundary of the model and matrix L chooses them from whole set of h M (h B = L T h M). It is important to note, that matrix S (e) chooses the field samples placed at the boundary of the outside region. In such a case, relation (12) represents electromagnetic properties of the selected region in a compact form. Note that matrix S ( e ) may be regarded as an equivalent of generalized admittance matrix of a region, except that γ rather than ω is a sweep parameter.We have arrived at the formulation that is analogous to the one considered in [6] in the context of MOR of RLC circuits. Namely, the ENOR algorithm was derived in [6] by considering a system of the form
where C, G and T are symmetric semi-positive matrices describing circuit reactive or resistive components, X is a vector of nodal voltages and J is a vector of current sources at ports whose location is described by the incidence matrix B and s is the Laplace variable s = jω (a complex angular frequency). It has been show that ENOR algorithm produces an orthogonal projection basis V, which spans the solution subspace for the assumed frequency band (s). The projection where X̂(s) = V T X(s), significantly reduces a dimension of the problem (dimX̂(s) << dimX(s)). As a result, one looses the information about the internal state variables, however the relation between output and input ports is preserved. The dimension of X̂(s) depends on the frequency bandwidth and on the assumed accuracy. If q is an order of the reduction and p is a number of ports, then relation (13) is approximated with respect to frequency up to its q moments (with a center point s 0) [6]. As a result of the projection with basis V, dim X̂(s) is equal p × q. Since dimX̂(s) << dimX, the inversion in (12) is performed for a significantly smaller matrix.Since (13) is analogous to (11) the same model order reduction technique can be applied to the wave propagation problem by replacing s with a propagation coefficient γ. Let V be a projection matrix obtained from ENOR algorithm for relation (11), then it can be reduced to the following form
where ĥ M = V T h M.Rather than forming the matrix valued transfer function for the modeled region, it is more convenient to perform the projection directly in curl operators (10) (see [9])
In this case, the problem has the same form as (10), but the number of variables is significantly smaller. Also in this case (like in the case concerning frequency) the accuracy of the approximation is high for a fixed range of propagation coefficient γ. To improve the accuracy or widen the bandwidth of γ a higher reduction order q is required.So far we have defined the input and output ports of a macromodel based on the field samples at the grid points adjacent to the boundary between the modeled region and a regular grid. When the mesh is fine the numbers of ports is large which has an adverse effect on the efficiency of MOR. This problem can be circumvented by applying discrete function expansion [1] to the field samples loacted at the boundary, prior to defining the ports of the macromodel.
In DFE the field samples in Yee’s mesh are replaced by several amplitudes of its harmonics. This approach can be treated as a second projection with the orthogonal basis [15]. In such a case, the field samples at the boundary are projected into smaller space of amplitudes (denoted by tilde)
where Q e and Q h are orthonormal matrices, whose columns corresponding to the boundary contain the discretized modal functions (in Yee’s mesh). As a result, Maxwell’s grid equations are projected two times - the first time using orthogonal basis (17) and the second time using basis obtained via ENOR algorithm. Due to DFE, the number of ports is considerably reduced: from the number of field samples at the boundary of the macromodel p to the number of expansion functions N.Finally, having performed all the projections we can eliminate the electric field arriving at the standard eigenvalue problem with γ 2 as an eigenvalue
where tilde and hat symbols represents the DFE and ENOR projections, respectively.3. Numerical results
To show the efficiency of the algorithm proposed in this paper we have analyzed three microstructured optical fibers. In the numerical tests we used the boundary conditions for modeling of open space (with 20 basis functions) and the scheme for a refractive index averaging (for a cells with a curved interface between different media) proposed in [1]. All calculations were performed in the MATLAB environment on a standard PC (i3 3.07GHz).
The first structure used in our tests (Fig. 4a) is a PCF (Photonic Crystal Fiber) with six circular holes (air) arranged in a hexagonal setting. The radius of the holes is r = 2.5μm, the pitch length is Λ = 6.75μm and the refractive index of the background medium is nbg = 1.45. In the analysis the vacuum wavelength λ = 1.45μm is assumed. Due to the symmetry of the structure only a quarter of the cross section (circle radius R = 10μm) is analysed. Hence, the boundary consist of a PEC (Perfect Electric Conductor) and/or PMC (Perfect Magnetic Conductor) at the symmetry planes and the radiation condition at the arc boundary [1]- see Fig. 4b. The numerical domain is divided into 200 cells along ρ direction and 160 along ϕ. For homogeneous regions (rings/layers 1 – 80 and 190 – 200) the DFE method with 20 basis functions is used. The macromodel is defined in the inhomogeneous part of the domain (rings 78 – 192). An application of the DFE technique at the boundary of the macromodel reduces the number of ports p from p = 960 (6 rings of 160 samples) to p = 120 (6 rings of 20 amplitudes). The propagation coefficient corresponding to an effective refractive index n 0 = 1.44 is assumed as a center of the “γ band”. The results obtained for different orders q are collected in Table 1 and compared to the ones obtained from standard FD scheme and FD scheme with DFE. The combination of the MOR and DFE techniques can significantly reduce the problem size which results in shorter computation time (about 14 times).
Table 1. Effective refractive index neff obtained for PCF from Fig. 4a. Eigenproblem solution time - round brackets, problem size - square brackets.
A more detailed analysis of the influence of the reduction order q and the number of basis functions N on the results is shown in Fig. 5. As expected, the accuracy of the calculations depends on the field distribution - for higher order modes the higher number N is required. However, the reduction order q must be grater than 5 in both cases.
Fig. 5 Error of the effective refractive index as a function of the reduction order and the number of basis function (it is evaluated in reference to the result obtained from standard FD method).
To test the algorithm accuracy Group Velocity Dispersion (GVD) of the fundamental mode is calculated. In Fig. 6 the results obtained for q = 6 and q = 7 are compared with the ones from the standard FD method and a very good agreement is achieved (error does not exceed 0.5% for q = 7).
In order to show the flexibility of the method the structure was distorted, the diameters of the holes are different (see Fig. 7). Despite lower symmetry of such structure the problem can be efficiently reduced. The results obtained for different modes with q = 6 and N = 20 are collected in Table 2. All the parameters of the simulation (size of the reduced problem, eigenproblem solution time) are the same as for undistorted structure.
Table 2. Effective refractive index neff obtained for distorted PCF shown in Fig. 7.
The next considered microstructure is a fiber presented in Fig. 8a. The angular-shaped holes have the inner radius r 1 = 1μm, the outer radius r 2 = 2μm and the angular width of 108°. The refractive index of the background medium is nbg = 1.44402362. In this case, the analysis is carried out for vacuum wavelength λ = 1.55μm. Also a shape of the numerical domain is different for this example - only a half of the cross section can be analyzed (circle radius R = 2.25μm) - see Fig. 8b. The numerical domain is divided into 140 cells along ρ direction and 210 along ϕ. Also in this case the DFE method with 20 basis functions is used for the homogeneous regions (rings 1 – 45 and 120 – 140). The macromodel is defined in the inhomogeneous part of the domain (rings 42 – 123). A propagation coefficient corresponding to the effective refractive index n 0 = 1.28 is assumed as a center of the “γ band”. The results obtained for different orders q are collected in Table 3. The projections reduced the number of state variables from almost 60000 to 3000 which significantly shortens the time of eigenvalue problem solution (from 8s to 0.5s).
Table 3. Effective refractive index neff obtained for PCF from Fig. 8a. Eigenproblem solution time - round brackets, problem size - square brackets.
The geometry of the last PCF is presented in Fig. 9a. The radius of the holes and the core is r = 2μm and the pitch length is Λ = 5μm. The refractive index of the background medium is nbg = 1.42, core nc = 1.45. The vacuum wavelength λ = 1.5μm is assumed. The discretization along ρ is M = 200 and along ϕ is N = 160. Similarly to the first fiber, the analysis is limited to a quarter of a circle of radius R = 8μm - see Fig. 9b. The DFE method with 20 basis functions is used for homogeneous regions (rings 1 – 80 and 190 – 200) and the macromodel is defined in inhomogeneous part of the domain (rings 78 – 192). A propagation coefficient corresponding to an effective refractive index n 0 = 1.41 is assumed as a center of the “γ band”. The results obtained for different orders q are collected in table 4. Also for this PCF the results well agree with the ones obtained from the standard FD method. The problem was reduced from over 60000 to just 3900 variables and the analysis was 16 times shorter.
Table 4. Effective refractive index neff obtained for PCF from Fig. 9a. Eigenproblem solution time - round brackets, problem size - square brackets.
In all examples the accuracy of the analysis depends on the reduction order q, which has to be high enough (q > 5). Small discrepancy in the results can be observed for propagation coefficients which differs from the center of the “γ band”. However, it is a consequence of the assumed approximation and to improve the accuracy a higher reduction order must be used.
4. Conclusion
In this paper a MOR technique was adopted to reduce the size of a propagation problem formulated in FDFD method. A verification of the proposed algorithm was performed on three different types of microstructured optical fibers. In all presented cases the ENOR algorithm combined with DFE method results in significant reduction of the problem size (up to 20 times) and the simulation time (up to 16 times).
Acknowledgments
This work has been funded by Polish Ministry of Science and Higher Education under grant 5407/B/T02/2010/38.
References and links
1. P. Kowalczyk, M. Wiktor, and M. Mrozowski, “Efficient finite difference analysis of microstructured optical fibers,” Opt. Express13, 10349–10359 (2005), http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-25-10349.
2. P. Kowalczyk and M. Mrozowski, “A new conformal radiation boundary condition for high accuracy finite difference analysis of open waveguides,” Opt. Express15, 12605–12618 (2007), http://www.opticsinfobase.org/abstract.cfm?URI=oe-15-20-12605. [CrossRef] [PubMed]
3. B. Moore, “Principal component analysis in linear systems: Controllability, observability, and model reduction,” IEEE Trans. Automat. Contr. , AC-26, 1732 (1981).
4. P. Feldmann and R. W. Freund, “Efficient Linear Circuit Analysis by Pade Approximation via the Lanczos Process,” IEEE Transactions on Computer-Aided Design , 14, 639–649 (1995). [CrossRef]
5. A. Odabasioglu, M. Celk, and L. T. Pileggi, “PRIMA: Passive Reduced-order Interconnect Macromodeling Algorithm,” 34th DAC, 58–65 (1997).
6. B. N. Sheehan, “ENOR: Model Order Reduction of RLC Circuits Using Nodal Equations for Efficient Factorization,” in Proc. IEEE 36th Design Automat. Conf., 17–21 (1999).
7. L. Kulas and M. Mrozowski, “Macromodels in the frequency domain analysis of microwave resonators,”Microwave and Wireless Components Letters, IEEE, 14, 94–96 (2004), http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1278378&isnumber=28582 [CrossRef]
8. L. Kulas and M. Mrozowski, “Multilevel model order reduction,” Microwave and Wireless Components Letters, IEEE, 14, 165–167 (2004), http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1291452&isnumber=28762 [CrossRef]
9. L. Kulas and M. Mrozowski, “A fast high-resolution 3-D finite-difference time-domain scheme with macromodels,” IEEE Trans. Microwave Theory Tech. 52, 2330– 2335 (2004). [CrossRef]
10. J. Podwalski, L. Kulas, P. Sypek, and M. Mrozowski, “Analysis of a High-Quality Photonic Crystal Resonator,” Microwaves, Radar & Wireless Communications. MIKON 2006. International Conference on, 793–796 (2006) http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4345301&isnumber=4345079.
11. A. C. Cangellaris, M. Celik, S. Pasha, and Z. Li,“Electromagnetic model order reduction for system-level modeling,” Microwave Theory Techniques, IEEE Transactions on, 47, 840–850 (1999), http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=769317&isnumber=16668. [CrossRef]
12. Y. Zhu and A. C. Canellaris, Multigrid Finite Element Methods for Electromagnetic Field Modeling (John Wiley & Sons, Inc., 2006). [CrossRef]
13. Z. Zhu and T. G. Brown, “Full-vectorial finite-difference analysis of microstructured optical fibers,” Opt. Express10, 853–864 (2002), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-17-853. [PubMed]
14. S. Guo, F. Wu, S. Albin, H. Tai, and R. S. Rogowski, “Loss and dispersion analysis of microstructured fibers by finite-difference method,” Opt. Express12, 3341–3352 (2004), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-15-3341.
15. L. Kulas, P. Kowalczyk, and M. Mrozowski, “A Novel Modal Technique for Time and Frequency Domain Analysis of Waveguide Components,” Microwave and Wireless Components Letters, IEEE, 21, 7–9 (2011), http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=5659497&isnumber=5680668. [CrossRef]
