Interior Penalty Discontinuous Galerkin Finite Element ... - IEEE Xplore

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 58, NO. 12, DECEMBER 2010

IV. CONCLUSION In this communication, a convenient form of a soft diffraction coefficient for a cylinder-tipped wedge was derived in order for it to be implemented into generic and multi-purpose ray-tracing algorithms. The diffraction coefficient is a UTD-based expression augmented by a rapidly converging series representing the Correction Field term. The numerical results have shown a high degree of accuracy for different cases of cylinder-tipped wedges. ACKNOWLEDGMENT A. C. Polycarpou would like to thank the members of COST action IC0603 for constructive comments and discussions on the topic presented in this communication.

4085

Interior Penalty Discontinuous Galerkin Finite Element Method for the Time-Dependent First Order Maxwell’s Equations Stylianos Dosopoulos and Jin-Fa Lee

Abstract—An interior penalty discontinuous Galerkin method is described with a conformal perfectly matched layer (PML) for solving the two first-order Maxwell’s equations in the time domain. Both central and upwind fluxes are studied in this work. In both cases, the proposed method is explicit and conditionally stable. Additionally, a local time-stepping strategy is applied to increase efficiency and reduce the computational time. Finally, numerical examples are presented to validate the method. Index Terms—Conformal PML, discontinuous Galerkin (DG), local time-stepping, Maxwell’s equations, time-domain (TD).

REFERENCES

I. INTRODUCTION

[1] T. K. Sarkar, Z. Ji, K. Kim, A. Medour, and M. Salazar-Palmar, “A survey of various propagation models for mobile communications,” IEEE Antennas Propag. Mag., vol. 45, no. 3, pp. 51–82, Jun. 2003. [2] J. Baldauf, S.-W. Lee, L. Lin, S.-K. Jeng, S. M. Scarborough, and C. L. Yu, “High frequency scattering from trihedral corner reflectors and other benchmark targets: SBR versus experiment,” IEEE Trans. Antennas Propag., vol. AP-39, no. 9, pp. 1345–1351, 1991. [3] R. G. Kouyoumjian and P. H. Pathak, “A uniform geometrical theory of diffraction for an edge in a perfectly conducting surface,” Proc. IEEE, vol. 62, no. 11, pp. 1448–1461, Nov. 1974. [4] J. B. Keller, “How dark is the shadow of a round ended screen?,” J. Appl. Phys., vol. 30, pp. 1452–1454, Sept. 1959. [5] R. G. Kouyoumjian and W. D. Burnside, “The diffraction of a cylindertipped half-plane,” IEEE Trans. Antennas Propag., vol. AP-18, no. 3, pp. 424–426, 1970. [6] W. Hallidy, “On uniform asymptotic Green’s functions for the perfectly conducting cylinder tipped wedge,” IEEE Trans. Antennas Propag., vol. AP-33, no. 9, pp. 1020–1025, 1985. [7] J. W. Yu and N. H. Myung, “TM scattering by a wedge with concaved edge,” IEEE Trans. Antennas Propag., vol. AP-45, no. 8, pp. 1315–1316, 1997. [8] J. J. Kim, H. J. Eom, and K. C. Hwang, “Electromagnetic scattering from a slotted conducting wedge,” IEEE Trans. Antennas Propag., vol. AP-58, no. 1, pp. 222–226, 2010. [9] K. C. Hwang, “Scattering from a grooved conducting wedge,” IEEE Trans. Antennas Propag., vol. AP-57, no. 8, pp. 2498–2500, 2009. [10] R. A. Ross and M. A. K. Hamid, “Scattering by a wedge with rounded edge,” IEEE Trans. Antennas Propag., vol. AP-19, no. 4, pp. 507–516, 1971. [11] G. Koutitas and C. Tzaras, “Rounded wedge diffraction coefficient including the effect of curvature discontinuity,” in Int. Symp. on Antennas and Propag., 2007. [12] Y. T. Chu, V. E. Anderson, H. H. Hubbel, Jr., and T. L. Ferrell, “Asymptotic solution for the diffraction of an electromagnetic plane wave by a cylinder-tipped half plane,” J. Opt. Soc. Am., vol. 73, no. 6, pp. 768–775, 1983. [13] M. A. Lyalinov, “Generalized Sommerfeld integral and diffraction in an angle-shaped domain with a radial perturbation,” J. Phys. A: Math. Gen., vol. 38, pp. L707–L714, 2005. [14] S. D. Weiner and S. L. Borison, “Radar scattering from blunted cone tips,” IEEE Trans. Antennas Propag., vol. AP-14, no. 6, pp. 774–781, Nov. 1971. [15] S. N. Karp, Diffraction by a Tipped Wedge New York Univ., Inst. Math. Sci., Div. EM Res. Rep. EM-52. [16] G. Tyras, Radiation and Propagation of Electromagnetic Waves. New York: Academic, 1969. [17] C. A. Balanis, Advanced Engineering Electromagnetics. New York: Wiley, 1989. [18] A. C. Polycarpou, Introduction to the Finite Element Method in Electromagnetics. London, U.K.: Morgan & Claypool, 2006.

The high level of complexity in real-life electromagnetic applications requires the development of sophisticated methods to solve the time-dependent Maxwell’s equations. Discontinuous Galerkin finite element methods are highly suitable for time-domain simulations. They support various types and shapes of elements, non-conformal meshes and non-uniform degrees of approximation. Furthermore, the resulting mass matrix is a block diagonal matrix and the method can lead into a fully explicit time-marching schemes. Finally, information exchange is required only between neighboring elements regardless of order and shape, which then leads to high efficiency in parallelization. In electromagnetism, DG methods were recently applied for the solution of the time dependent Maxwell’s equations [1]–[4]. The main contribution of this work is the combination of a conformal perfectly matched layer with local time-stepping techniques in DGTD [3] and the numerical study of the convergence of DGTD methods based on centered and upwind fluxes for the first order Maxwell’s equations. The remainder of the communication is organized as follows. In Sections II, III we present the space and time discretization of the interior penalty (IP) [5], [6] DG formulation for the first order Maxwell’s equations in the time-domain. Next, in Section IV, we present how to incorporate a conformal perfectly matched layer in DGTD and discuss the issues concerning the stability and performance of the perfectly matched layer (PML). In Section V we numerically study the convergence rate for both central and upwind fluxes and compare with the theoretical estimates derived in [1], [2] accordingly. Finally, in Section VI. the validation of the method is presented through numerical examples and the local time-stepping strategy is employed in multi-scale applications to illustrate the potential and the advantages of the strategy.

II. DGTD FORMULATION-SPACE DISCRETIZATION Following the interior penalty approach described in [5], [6] we can derive the DGTD formulation in (1). Assuming constant material Manuscript received August 03, 2009; revised May 03, 2010; accepted May 18, 2010. Date of publication September 23, 2010; date of current version November 30, 2010. The authors are with the Electroscience Lab, Departement of Electrical and Computer Engineering, The Ohio State University, Columbus, OH 43212 USA (e-mail: [email protected]; [email protected]). Color versions of one or more of the figures in this communication are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TAP.2010.2078445

0018-926X/$26.00 © 2010 IEEE

4086


(Phii)nm = 0

properties within each element we can define the following finite-di3 mensional discrete trial space: Vhk = fv 2 [L2 ( )] : vjK 2 3 k [P (K )] ; 8K 2 Th g. The final formulation can be formally stated as

F ind (H; E)

2 Vhk 2 Vhk

such that 1 @ H w1 r2E+ d

@t

1 @ E d + ffvgg 1 [ H] ds 0 v1 r2H0

@t

F 0

F

(Peij )nm =

ffwgg 1 [ E] ds + e

F

(Phij )nm =

[ v] 1 [ E] ds

+ f [ w] 1 [ H] ds = 0: F 8(w; v) 2 Vhk 2 Vhk :

(1)

We define the tangential trace and projection (“components trace”) op^ i 2 ui j@ and erators, (1) and (1) respectively, as (ui ) = n (ui ) = n^ i 2 (ui 2 n^ i )j@ , and n^ i is the boundary normal pointing out of the element Ki . Moreover, ffugg = ( (ui ) + (uj ))=2, [ u] = (ui ) + (uj ) and [ u] = (ui ) 0 (uj ). For the choice e = f = 0 one obtains a conservative formulation, but with a suboptimal (O(hp )) rate of convergence as discussed in [2]. If we choose, e = 1=(2Z0 ) and f = 1=(2Y0 ) with Z0 = (1=2)( i =i + j =j ), Y0 = (1=2)( i =i + j =j ), we get a slightly lossy formulation. However, in this case an optimal (O(hp+1 )) rate of convergence is obtained as proven in [1]. III. DGTD FORMULATION-TIME DISCRETIZATION Let us expand the electric and magnetic fields within element Ki in terms of basis functions w; v 2 Vhk as E(hi) (r; t) = E(r; t)jK dn ein (t)vin (r) and Hh(i) (r; t) = H(r; t)jK d n hin (t)win (r) where di are the degrees of freedom in element Ki . Separating the w, v testing in (1) we get in matrix form a semi-discrete system within each element

@ ei =Se hi 0 Feii hi + ePeiiei 0 Feij hj + ePeij ej @t @h M i = 0 Sh ei + Fhii ei + f Phii hi + Fhij ej + f Phij hj @t

M

(M )nm = (Se )nm = (Sh )nm =

K K K K

vin 1 (i 1 vim )dK win 1 (i 1 wim )dK vin 1 (r 2 wim )dK win 1 (r 2 vim )dK

(Feii)nm = 12 F 1 (Fhii)nm = 2 F 1 ij (Fe )nm = 2 F 1 (Fhij )nm = 2 F ii (Pe )nm = 0

F

(vin ) 1 (wim ) (win ) 1 (vim )ds (vin ) 1 (wjm )ds (win ) 1 (vjm )ds (vin ) 1 (vim )ds

F

(win ) 1 (wim )ds

(vin ) 1 (vjm )ds (win ) 1 (wjm )ds

are the local matrices. The system of ODEs (2), (3) is discretized in time with a leap-frog scheme, which is second order accurate. The electric field unknowns are evaluated at tn = nt and the magnetic field unknowns are evaluated at tn+(1=2) = (n + (1=2))t. The first derivatives will be approximated using the central difference method. Moreover, for the two extra penalty terms arising from the upwind flux formulation, we use a backward approxin+(1=2) n+(1=2) mation ei(j ) eni(j ) and hni(+1 . If an average j ) hi(j ) n+(1=2) approximation was used i.e. (ei(j ) (ein(j ) + eni(+1 j ) )=2 and n+(1=2) n+(3=2) n +1 hi(j ) (hi(j ) + hi(j) )=2), then the system will become globally implicit due to the coupling terms from the neighboring elements. Therefore, the backward approximations in the upwind flux formulation are necessary if we want the time marching scheme to remain explicit. In this way, the fully discretized local system of equations is written as

M eni +1 = M + etPeii eni + t Se 0 Feii hi

n+

+ etPeij enj = M + ftPhii hni + + t 0Sh + Fhii + tFijh enj +1 + ftPhij hnj + : n+

n+ M hi

0 tFij e hj

(4)

eni +1 (5)

Let us note that, the leap frog scheme is second order accurate. The overall error can be written as: e = O(hp+1 )+ O(t2 ) O(hp+1 )+ O(h2 =c). Therefore to maintain the optimal O(hp+1 ) rate for higher order basis in space a higher order time discretization like LF-4 is also needed.

(2)

A. Stability Analysis and Local Time-Stepping (LTS)

(3)

The resulting update scheme is conditionally stable. A stability condition ispderived in [2], [3]. For first-order DGTD: 8i; 8j 2 N (i), ci ti [(4 5=3) + (8=3)max( i =j ; i =j )] 0, where h is the mesh size. The problem under consideration is the one of a TEM wave propagating in a parallel plate waveguide. In this case k = (0kx ; 0; 0), E = (0; 0;Ez ) and H = (0;Hy ; 0). The excitation is given by Einc = (0; 0; 1)ej!t jk r and Hinc = (1=)(0; 1; 0)ej!t jk r . The material properties are r = 1 and r = 1. The waveguide was terminated with a first order Silver-Muller absorbing boundary condition ^ 2 E = 0cn^ 2 n^ 2 H and n^ 2 H = cn^ 2 n^ 2 E), which is (i.e. n 0

1

0

1

exact for this problem. In Fig. 4, we present numerical results for the convergence studies using central and upwind fluxes for both electric and magnetic fields. The results show that central flux provides a suboptimal convergence a ( p ) convergence rate for central flux and a ( p+1 ) convergence rate for an upwind flux, as predicted by the theoretical estimate derived in [1], [2] accordingly.

Oh

Oh


4089

= 30m

Fig. 6. RCS (dB) of a coated sphere with inner radius a : and outer : radius b at 300 MHz. (Solid line-DGTD, dashed line-Mie series). (a) H-plane and (b) E-plane.

= 3 25 m

VI. NUMERICAL RESULTS In all of the examples we used tetrahedral meshes and first order polynomials for the approximation of the fields inside each element. Therefore, we have di = 12 degrees of freedom inside each element. For all numerical examples, an upwind flux formulation was used.

TABLE I ELEMENT PARTITIONING BY CLASSES FOR THE F16-AIRCRAFT

A. Coated Sphere In this example, we calculate the RCS of a coated sphere and compare it with the exact solution. We simulated a coated sphere with inner radius of a = 3:0 m, the outer radius is b = 3:25 m and the coating has r = 2:0 and r = 1:0. Due to the limited computational resources (32 GB of RAM), the generated mesh has an average edge length lavg (=7) at 300 MHz in free space. The PML/air boundary is set 0.5 m from the coating. The PML has a constant profile of = 0:02 and a thickness of 0.4 m. The mesh consisted of 2,242,037 elements and required 26.5 GB of memory. The coated sphere is illu^ minated by a Gaussian pulse Einc = E0 e0[t0t 0k1(r0r )=c] = with E0 = (1; 0; 0), k^ = (0; 0; 1), t0 = 5:0 ns, = 0:6 ns. By employing the Fourier transform of the time-dependent electric and magnetic currents on the sphere’s surface, we can evaluate the RCS at multiple frequencies. Fig. 6 shows the computed RCS at the frequency of 300 MHz. A good agreement is observed between the numerical solution and the analytical solution provided by the Mie series. The RCS error defined by (15) was computed to be RCSerror = 9:02 2 1002

RCSerror =

Fig. 7. (a) Surface mesh of the F16-aircraft. (b) Geometry of F16-aircraft together with the PML truncation. The conformal PML as well as the glass canopy and dielectric radome are shown accordingly.

2 0 0

j (; )0 (; )j2 sindd : 2 2 j ( ; )j sindd 0 0 M ie

N um

M ie

(15)

B. Scattering From a F-16 Aircraft In this example, we consider the numerical simulation of the scattering from an F-16 aircraft. The aircraft is illuminated by a Gaussian ^ ^= pulse Einc = E0 e0[t0t 0k1(r0r )=c] = with E0 = (0; 0; 1), k (1; 0; 0), t0 = 4:0 ns, = 0:6 ns and r0 = (8:78; 0; 0) m. A partial view of the surface mesh is shown in Fig. 7. Let us note, that the modeling of the F-16 also includes the dielectric radome and the glass canopy as shown in Fig. 7. The generated mesh has an average edge length (=5) at 300 MHz, which is the 3-dB frequency of the incident pulse. The total number of tetrahedra is 1,656,676 and the total number of DOF is 52,641,744 with 13,213,392 DOF in the PML region. The PML is set to be 4-layers thick and has a constant profile of = 0:02. In the outer surface of the PML region, a first order ABC was applied to reduce the reflection due to the finite thickness of the PML. The mesh generated in this case is an unstructured locally refined

TABLE II CPU TIME GAIN OBTAINED WITH THE LOCAL TIME-STEPPING VERSUS THE LEAP-FROG SCHEME

mesh with strong element-size disparities, which leads to a minimum time step of tmin = 7:50 2 10013 s and a maximum time step of tmax = 5:33 2 10011 s. The application of the local time-stepping strategy is shown in Table I. The elements are partitioned into 4 classes. Only 0.015% of the elements are the really small elements, that would have dominated in the stability condition and increased significantly the overall CPU time. Moreover, a simple parallelization for shared memory architecture using Open MP parallelization was implemented. This simulation was performed on a two quad-core Xeon processors with 32 GB RAM machine, and the number of processors was Np = 6. For this example the application of the local time-stepping will provide a gain of approximately 15 times in the CPU time as shown in Table II. Let us note, that in the non local time stepping case the CPU time is estimated time based on the average iteration time. From Table II, one can clearly see the advantages of the local time-stepping algorithm for meshes with strong element-size disparities, which are commonly found in real-life applications. Finally, the field distributions for the electric field at different times is given at Figs. 8, 9. One can notice the significant scattering created by the air intake as well as by the nose and the cockpit regions.

4090


Fig. 8. Electric field distribution on the surface of the aircraft at different time locations. (a) Top view. (b) Bottom view.

[6] P. Houston, I. Perugia, A. Schneebeli, and D. Schotzau, “Interior penalty method for the indefinite time-harmonic maxwell equations,” Numer. Math., vol. 100, no. 3, pp. 485–518, 2005. [7] J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys., vol. 114, no. 2, pp. 185–200, 1994. [8] Z. Sacks, D. Kingsland, R. Lee, and J.-F. Lee, “A perfectly matched anisotropic absorber for use as an absorbing boundary condition,” IEEE Trans. Antennas Propag., vol. 43, no. 12, pp. 1460–1463, Dec. 1995. [9] S. Gedney, “An anisotropic perfectly matched layer-absorbing medium for the truncation of fdtd lattices,” IEEE Trans. Antennas Propag., vol. 44, no. 12, pp. 1630–1639, Dec. 1996. [10] F. L. Teixeira and W. C. Chew, “Analytical derivation of a conformal perfectly matched absorber for electromagnetic waves,” Microwave Opt. Technol. Lett, vol. 17, pp. 231–236, 1997. [11] S. Abarbanel, D. Gottlieb, and J. S. Hesthaven, “Long time behavior of the perfectly matched layer equations in computational electromagnetics,” J. Sci. Comput., vol. 17, no. 1–4, pp. 405–422, 2002. [12] J. Niegemann, M. König, K. Stannigel, and K. Busch, “Higher-order time-domain methods for the analysis of nano-photonic systems,” Photon. Nanostructures—Fundamentals Applicat., vol. 7, no. 1, pp. 2–11, 2009.

A Domain Decomposition Technique for Efficient Iterative Solution of Nonlinear Electromagnetic Problems Giacomo Guarnieri, Giuseppe Pelosi, Lorenzo Rossi, and Stefano Selleri

Fig. 9. Snapshots of the magnitude of the electric field distribution in the computational domain. The F-16 aircraft is illuminated by a Gaussian pulse on-nose incidence.

VII. CONCLUSION We have presented an IP-DGTD with a conformal PML for the domain truncation. The numerical convergence of both central and upwind flux formulations were studied. Moreover, the late time instability and performance of the PML were studied by means of numerical experiments. Finally, a local time-stepping technique was employed to reduce the total CPU time and improve the efficiency in multi-scale applications. ACKNOWLEDGMENT The authors thank Ansoft Corporation for their support of this work.

REFERENCES [1] J. S. Hesthaven and T. Warburton, “Nodal high-order methods on unstructured grids,” J. Comput. Phys., vol. 181, no. 1, pp. 186–221, 2002. [2] L. Fezoui, S. Lanteri, S. Lohrengel, and S. Piperno, “Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous maxwell equations on unstructured meshes,” ESAIM: M2AN, vol. 39, no. 6, pp. 1149–1176, 2005. [3] E. Montseny, S. Pernet, X. Ferriéres, and G. Cohen, “Dissipative terms and local time-stepping improvements in a spatial high order discontinuous Galerkin scheme for the time-domain Maxwell’s equations,” J. Comput. Phys., vol. 227, no. 14, 2008. [4] S. Gedney, C. Luo, J. Roden, R. Crawford, B. Guernsey, J. Miller, and E. Lucas, “A discontinuous Galerkin finite element time domain method with PML,” in Proc. IEEE AP-S Int. Symp., Jul. 2008, pp. 1–4. [5] D. N. Arnold, “An interior penalty finite element method with discontinuous elements,” SIAM J. Numer. Analy., vol. 19, no. 4, pp. 742–760, 1982.

Abstract—An innovative method for fast analysis of electromagnetic problems comprising nonlinear dielectrics is proposed. The method is based on the combination of finite elements and domain decomposition methods. This latter technique allows for the separation of the problem in multiple subproblems which can be solved separately. By appropriately defining one of such subdomains as containing all the nonlinear dielectrics it is possible to restrict the application of an iterative algorithm for the numerical solution of nonlinear equations only to this latter, smaller, domain. Numerical results showing the accuracy and efficiency of this technique are presented in few 2D cases. Index Terms—Domain decomposition, finite element method, nonlinear media, nonlinear wave propagation.

I. INTRODUCTION The presence of nonlinear media usually makes the computation of the solution of Maxwell equations very demanding, yet there are many practical electromagnetic engineering problems in which nonlinear materials are essential, hence the scientific community has deeply investigated this problem, and a brief summary of these efforts is reported here. Few exact analytical results are available, for example for the scattering of a TM plane wave from a nonlinear thin film and for few other cases [1]. Approximate solutions, both analytical and numerical have wider applications. Among the earliest approximate techniques applied to nonlinear electromagnetics there is the perturbation technique used Manuscript received November 20, 2009; revised April 07, 2010; accepted May 14, 2010. Date of publication September 23, 2010; date of current version November 30, 2010. G. Guarnieri is with the BU Radar Systems, Dep. of Antennas, Selex Galileo S.p.A., I-50013 Campi Bisenzio, Florence, Italy. G. Pelosi, L. Rossi and S. Selleri are with the Department of Electronics and Telecommunications, University of Florence, I-50134 Florence, Italy (e-mail: [email protected]). Digital Object Identifier 10.1109/TAP.2010.2078462

0018-926X/$26.00 © 2010 IEEE