CC BY
32
UDC 519.63
Parallel Second Order Finite Volume Scheme for Maxwell's Equations with Discontinuous Dielectric Permittivity on
Structured Meshes
T. Z. Ismagilov
Department of Information Technology Novosibirsk State University 2, Pirogova str., Novosibirsk, Russia, 630090
A second order finite volume scheme on structured meshes is presented for numerical solution of time dependent Maxwell's equations with discontinuous dielectric permittivity. The scheme is based on approaches of Godunov, Van Leer and Lax Wendroff and employs a special technique for gradient calculation near dielectric permittivity discontinuities. The scheme was tested for problems with linear and curvilinear discontinuities. Test results demonstrate second order of convergence and support second order of approximation in space and time. A parallel implementation of the scheme based on geometric decomposition was developed. Computational region was partitioned into subregions. Computations in each subregion were carried out independently using halo cells. Test results indicate linear scalability. Parallel implementation was applied to modelling photonic crystal devices. Computational results for photonic crystal waveguide with a bend correctly confirm bend configurations and frequencies with zero reflection.
Key words and phrases: Maxwell's equations, Godunov scheme, finite volume, discontinuous permittivity, second order, photonic crystals, waveguides.
1. Introduction
Finite difference time domain method based on structured cartesian grids is arguably the most popular method for numerical solution of Maxwell's equations [1]. It is second order accurate in space and time for media with constant dielectric permittivity but has reduced order of approximation for media with dielectric permittivity discontinuities. Recently several finite volume schemes on unstructured meshes were suggested that are second order accurate in space and time even for media with dielectric permittivity discontinuities [2,3].
For many problems the use of unstructured meshes is not necessary. Structured meshes offer several advantages. They can be generated using trivial algebraic algorithms and schemes on structured meshes can be easily parallelized.
In this paper we suggest a second order finite volume scheme on structured meshes for numerical solution of Maxwell's equations with discontinuous dielectric permittivity with parallel implementation. The scheme uses Godunov flux approximation [4] and approaches of Van Leer [5] and Lax-Wendroff [6] to increase order of approximation. The key idea of the scheme is to use stencils for gradient approximation that don't cross dielectric permittivity discontinuity. Scheme was tested for linear as well as curvilinear discontinuities. Calculation results confirm second order of approximation. Parallel implementation was developed using OpenMP. Test results indicate linear scalability. Scheme was applied to modelling photonic crystal waveguides [7,8].
2. Maxwell's Equations
The system of two-dimensional Maxwell's equations for TM case can be written in vector form as
d d d
dtU + dx1 1 + dx2 2 0'
Received 30th December, 2013.
Ismagilov T. Z. Parallel Second Order Finite Volume Scheme for Maxwell's .
221
where U = (D3,Bi, B2) is conservative variables vector, Fi = (-H2,0, -E3) and F2 = (Hi,E3, 0)4 are flux vectors.
In the above formulas E is electric field, H — magnetic field, D = eE — electric induction, B = ^H — magnetic induction, e — dielectric permittivity, ^ — magnetic permeability. In this paper we assume ^ = 1. The system of Maxwell's equations can also be written using flux variables V = (E3, Hi,H2)t related to conservative variables by U = QV as
d , d , d
where
Q-V + Ai— V + A^— V = 0, at oxi ox2
(2)
e 0 0 Q = | 0 » 0 0 0 p
0 0 -1 Ai = | 0 0 0 -10 0
¿2 =
010 1 0 0 000
(3)
3. Numerical Scheme
By integrating the system of Maxwell's equations over a quadrilateral cell Ci with edges r k assuming constant dielectric permittivity in the cell an integral conservation law can be obtained
Q^ f VdQ + ^ f (mFi + n2F2) dr = 0, (4)
Ci r k
where (ni,n2) is a unit normal. For approximation of this integral conservation law consider a finite volume Godunov scheme
Vn+i - V " A , , QQC% V-^ + ^ s*Ci Fk = 0, (5)
T k=i
where QCi is volume of the ith cell, sC — length of its kth edge, t — time step. Flux F is calculated using exact solution to the Riemann problem F = l (Xr) + A-Vr (Xr) where Vl, R(Xr) are interpolations of V from two neighboring cells L and R on edge center Xr and matrices A+ and A- can be written as
1
±^elSR —eRn 2 -—eRm
A± = —=-— 1 -ein2 ±«2 Tnin2 I . (6)
^ + ^V -VLLm Tnin2 ±n2 '
The scheme will have second order of approximation in space and time if values at the edge centers are calculated with second order of approximation. Such values can be obtained using interpolation
V l, R(Xr) = V(Xl,r) + g(X L,R)(Xr - X l,r) - 2Q-i £ AiIV(X l,r), (7)
3=i 3
where X l,r are barycenters of neighboring cells, as long as the derivatives are approximated with first order [5,6].
Derivatives of V are approximated using the least squares method. Stencil for derivative approximation consists of a cell and 8 adjacent cells. If a cell is next to
dielectric permittivity discontinuity the stencil is shifted one cell away from discontinuity. So that all the cells in the stencil have the same dielectric permittivity as in the cell for which derivatives are calculated. Derivatives of V are limited by interpolating cell values on adjacent vertices, taking arithmetic average of interpolated values as vertex value and finally calculating limited gradients in a cell by applying least squares method to cell vertex values.
Parallel implementation using geometric decomposition requires data from up to 4 halo cells outside of subregion to find fields at a new time step in subregion. Parallel implementation was programmed using OpenMP. Computational region was partitioned in vertical subregions. Computation of new values in a subregion was programmed as a separate method. This method in a cycle row by row reads data from computational region and halo cells and computes the next row of electromagnetic fields at a new time step. In this way only the data from the rows necessary to update next row is stored in memory and not data from the whole subregion and all halo cells.
Figure 1. Computational mesh and electric field distribution
4. Numerical Experiments
To verify order of approximation of the scheme we considered problem of interaction of plane electromagnetic wave with a dielectric cylinder. The cylinder dielectric material had dielectric permittivity 2. A sequence of five meshes was used. Structured 40 by 40 mesh consisting of 1600 quadrilaterals and contour plot showing distribution of electric field at time 2.0 obtained using 80 by 80 mesh consisting of 6400 quadrilaterals are shown on Fig. 1. Maximum errors for different meshes are presented in Table 1. Error behavior indicates second order of convergence and supports second of approximation of the suggested scheme. Parallel implementation showed linear scalability. For example simulation using 640 by 640 mesh took 1009 and 507 seconds using 2 and 4 threads respectively.
To demonstrate scheme potential for solving real problems we considered problem of pulse propagation inside photonic crystal waveguide. Waveguide had a square lattice with cylindrical elements. All the parameters were from [8]. Reflection spectra obtained for two waveguide configurations are shown on Fig. 2 and show presence of one frequency with zero reflection for one configuration.
Ismagilov T. Z. Parallel Second Order Finite Volume Scheme for Maxwell's . . . 223
0.2 0.15
0.1
Figure 2. Reflection coefficients for photonic crystal waveguides
Table 1
Max. error cylinder
N Max. 62 Convergence order
40 0.051409 —
80 0.013267 1.95
160 0.003348 1.99
320 0.000840 2.00
640 0.000210 2.00
5. Conclusion
We have presented and tested a second order finite volume scheme on structured meshes for Maxwell's equations with discontinuous dielectric permittivity. Computational results for test problems confirm second order of approximation of the proposed scheme. Spectra obtained for photonic crystal waveguides agree well with the results obtained by other theoretical and computational approaches [7,8]. Parallel implementation of the scheme indicates linear scalability.
References
1. Yee K. S. Numerical Solution of Initial Boundary Value Problems Involving Maxwell's Equations in Isotropic Media // IEEE Trans. Antennas Propagat. — 1966. — Vol. 14. — Pp. 585-589.
2. Ismagilov T. Z. A Second-Order Scheme for Maxwell's Equations with Dielectric Permittivity Discontinuities and Total Field-Scattered Field Boundaries // Int. J. Comput. Math. — 2012. — Vol. 89, No 10. — Pp. 1378-1387.
3. Ismagilov T. Z. Second Order Scheme for Maxwell's Equations with Discontinuous Electromagnetic Properties // Lect. Notes Comp. Sci. — 2012. — Vol. 7125. — Pp. 227-233.
4. Годунов С. К. Разностный метод численного расчета разрывных решений уравнений гидродинамики // Математический сборник. — 1959. — Т. 47(89), № 3. — С. 271-306. [Godunov S.K. A Difference Method for Numerical Calculation of Discontinuous Solutions of the Equations of Hydrodynamics // Mat. Sb. (N.S.) — 1959. — Vol. 47(89), No3. — Pp. 271-306. — (in russian). ]
5. van Leer B. Towards the Ultimate Conservative Difference Scheme. V. A Second-Order Sequel to Godunov's Method // J. Comput. Phys. — 1979. — Vol. 32, No 1. — Pp. 101-136.
6. Lax P. D., Wendroff B. Systems of Conservation Laws // Commun. Pure Appl. Math. — 1960. — Vol. 13. — Pp. 217-237.
7. Ismagilov T. Z., Kuzmin A. I. Modelling Photonic Crystal Waveguides using Finite Volume Method // Days on Diffraction 30 May - 3 June 2011. St. Petersburg. — 2011. — Pp. 87-91.
8. A. Mekis, J. C. Chen, I. Kurland et al. // Phys. Rev. Lett. — 1996. — Vol. 77. — Pp. 3787-3790.
УДК 519.63
Параллельная конечно-объёмная схема второго порядка для уравнений Максвелла с разрывной диэлектрической проницаемостью на структурированных сетках
Т. З. Исмагилов
Факультет информационных технологий Новосибирский государственный университет ул. Пирогова, д. 2, Новосибирск, Россия, 630090
Предлагается схема второго порядка на структурированных сетках для численного решения нестационарных уравнений Максвелла с разрывной диэлектрической проницаемостью. Схема основана на подходах Годунова, Ван Леера и Лакса-Вендрофа и использует специальный подход к вычислению градиентов около разрывов диэлектрической проницаемости. Схема была проверена на задачах с линейными и криволинейными разрывами диэлектрической проницаемости. Результаты расчётов показывают второй порядок сходимости и подтверждают второй порядок аппроксимации схемы по времени и пространству. Используя метод геометрической декомпозиции, была разработана параллельная реализация схемы. Вычислительная область разбивалась на подобласти. Расчёты в каждой подобласти проводились независимо, используя дополнительные ячейки. Результаты расчётов подтверждают линейную масштабируемость. Параллельная реализация была применена для моделирования фотонно-кристаллических устройств. Результаты расчётов для фотонно-кристаллического волновода с изгибом правильно предсказывают конфигурации и частоты с нулевым отражением.
Ключевые слова: уравнения Максвелла, схема Годунова, метод конечных объёмов, разрывная диэлектрическая проницаемость, второй порядок, фотонные кристаллы, волноводы.
