Discontinuous Galerkin Time-Domain solution of ... - CiteSeerX

34 downloads 4209 Views 888KB Size Report
then be discretized similarly, through dot-multiplication by each local basis function ϕij ..... dipole located at the center of a [0, 1]3 domain in the free space.
Discontinuous Galerkin Time-Domain solution of Maxwell’s equations on locally-refined nonconforming Cartesian grids N. Canouet, France Telecom R&D, La Turbie, France, L. Fezoui and S. Piperno, INRIA-CERMICS, BP93, 06902 Sophia Antipolis Cedex, France

Keywords discontinuous Galerkin methods, centered fluxes, leap-frog time scheme, L 2 stability, energy methods, non-conforming grid Abstract The use of the prominent FDTD method for the time domain solution of electromagnetic wave propagation past devices with small geometrical details can require very fine grids and can lead to very important computational time and storage. Our purpose is to develop a numerical method able to handle possibly non-conforming locally refined grids, based on portions of Cartesian grids in order to use existing preand post-processing tools. We present a Discontinuous Galerkin method to solve the three-dimensional time-domain Maxwell’s equations on conforming or non-conforming orthogonal grids. The method is based on a centered mean approximation for surface integrals and a second-order leap-frog scheme for advancing in time. The dispersion error is partially analyzed in the cubic uniform mesh case. We show that the divergence of fields are weakly preserved on conforming grids. In the most general case, we prove that the resulting scheme is stable and that it conserves a discrete analog of the electromagnetic energy. For non-conforming grids, the local sets of basis functions are enriched at subgrid interfaces in order to get rid of possible spurious wave reflections. Numerical simulation on realistic configurations are promising: the numerical method proposed here makes it possible to handle devices with extremely small details. Further investigations are possible with different, higher-order discontinuous finite elements. Structured Abstract Research paper Purpose. The use of the prominent FDTD method for the time domain solution of electromagnetic wave propagation past devices with small geometrical details can require very fine grids and can lead to very important computational time and storage. Our purpose is to develop a numerical method able to handle possibly non-conforming locally refined grids, based on portions of Cartesian grids in order to use existing pre- and post-processing tools. Design/methodology/approach. We build a Discontinuous Galerkin method based on bricks and proof its stability, accuracy and efficiency. Findings. We find that it is possible to conserve exactly the electromagnteic energy and weakly preserves the divergence of the fields (on conforming grids). For non-conforming grids, the local sets of basis functions are enriched at subgrid interfaces in order to get rid of possible spurious wave reflections. Research limitations/implications. Although the dispersion analysis is incomplete, the numerical results are really encouraging and we show the proposed numerical method makes it possible to handle devices with extremely small details. Further investigations are possible with different, higher-order discontinuous finite elements. Originality/value. We claim this paper can be of great value for people wanting to migrate from FDTD methods to more uptodate time-domain methods, while conserving existing pre- and post-processing tools.

1

Introduction The modeling of systems involving electromagnetic waves is now widely done through the solution of the time-domain Maxwell equations on space grids. Although the Finite Difference Time-Domain (FDTD) methods based on Yee’s scheme (Yee 1966) are still prominent, many different types of methods have been proposed which are based on unstructured meshes and can deal with complex geometries, like Finite Element Time-Domain (FETD) methods (Lee & Sacks 1995, Yioultsis, Kantartzis, Antonopoulos & Tsiboukis 1998, Elmkies & Joly 1997, Joly & Poirier 2000), mimetic methods (Hyman & Shashkov 1999), Finite Volume Time-Domain (FVTD) methods (Holland, Cable & Wilson 1991, Cioni, Fezoui, Anne & Poupaud 1997, Remaki 2000, Lohrengel & Remaki 2002) and Discontinuous Galerkin Time-Domain methods, which enjoy a renewed favor nowadays and are now used in many and various applications (Remaki & Fezoui 1998, Cockburn, Karniadakis & Shu 2000) since they can achieve a high order of accuracy by simply choosing suitable basis functions, including spectral elements (Kopriva, Woodruff & Hussaini 2000), either on tetrahedral meshes using Lagrange polynomials (Hesthaven & Teng 2000, Hesthaven & Warburton 2002, Hesthaven & Warburton 2004) or on hexahedral meshes using products of Legendre polynomials. The existing software are mostly based on upwind or partially upwind fluxes, along with multi-step low-storage Runge-Kutta time-schemes (less often on a leap-frog scheme and centered fluxes, a non-diffusive combination which was first proposed by Remaki (Remaki 2000) to achieve long time wave propagation) and remain highly parallelizable. One can notice here that, although Runge-Kutta schemes are robust and stable, no stability proof is available for these kind of schemes on arbitrary meshes (stability is proved for FVTD methods on either Cartesian grdis (Remaki 2000) or unstructured meshes (Piperno, Remaki & Fezoui 2002)). At the same time, one of the most important properties which should be aimed at is the conservation of a discrete analog of the electromagnetic energy (if there is no electric conductivity). This cannot be obtained with Discontinuous Galerkin methods based on upwind fluxes (Cioni et al. 1997, Warburton 2000, Kopriva et al. 2000), although upwind fluxes lead to more robust codes, particularly for frequency-domain computations (low-frequency stabilization of centered fluxes have been proposed though (Hesthaven & Warburton 2004)). Another important property is the preservation of the divergence relations in the absence of sources: the electric and magnetic fields should remain solenoidal throughout the computation. Much work has recently been done in that direction, including divergencefree basis functions (Cockburn, Li & Shu 2004). Increasing numbers of numerical simulations of electromagnetic and electronic systems are being made on unstructured (mostly tetrahedral) grids rather than Cartesian grids and using DGTD methods rather FDTD methods. However, FDTD simulations and software are still prominent, maybe mainly because many preprocessing and post-processing software has been designed for Cartesian grids. It is clear that Yee’s scheme is very efficient for the propagation in a homogeneous medium (good phase properties, easy implementation) and as long as the geometry and the physical data (currents, material interfaces, etc..) remain simple. In presence of complex structures, one may still use the Yee’s scheme with very fine uniform grids but the solution becomes very costly in time and memory and sometimes impossible to compute. An alternative is to locally refine the grid, only where it is needed and then join the resulting subgrids in a conforming or non-conforming way. It is well known that the increasing of dispersion errors and the possible instability are the major drawbacks when applying explicit schemes to locally refined grids. Another possible source of inaccuracy are the possible spurious reflections appearing at interfaces between Cartesian subgrids. In this paper, we present a numerical method we have designed, which is aimed at (1) being efficient on Cartesian grids, (2) being able to handle locally refined grids, possibly in a non-conforming way, (3) conserving exactly a global electromagnetic energy, (4) being similar to recently developed approaches using Discontinuous Galerkin methods (Cockburn et al. 2000, Hesthaven & Warburton 2002, Hesthaven & Warburton 2004, Piperno & Fezoui 2003), including divergence-free approaches (Cockburn et al. 2004). It combines a centered flux approximation with a leap-frog time integration (resulting in a non diffusive scheme), the choice of a new set of basis functions well adapted to Cartesian grids, in particular concerning the exact preservation of divergence constraints. Also, we propose to increase the degree of the approximation at non-conforming interfaces. We will also show that the new scheme is stable, non diffusive (a discrete energy is conserved) and is able to work on highly refined grids without any

2

noticeable dispersion and reflection effects, including in the case of long time calculations.

The Discontinuous Galerkin Time-Domain (DGTD) method We consider the three-dimensional time-domain Maxwell’s equations for heterogeneous anisotropic linear media with source ~j (current density). The electric permittivity tensor ε¯(x) and the magnetic ¯(x) are symmetric positive definite, varying in space and uniformly bounded. The permeability tensor µ ~ = t(Ex , Ey , Ez ) and the magnetic field H ~ = t(Hx , Hy , Hz ) verify the following Maxwell’s electric field E equations, where divergence constraints have been dropped:  ~ ∂E   ~ H ~ − ~j,  ε¯ = curl ∂t (1) ~  ∂H   µ ~ ~ ¯ = −curl E. ∂t

These equations are set and solved on a bounded domain Ω of R3 . On the regular domain boundary ∂Ω ~ = ~0) (of unitary outwards normal ~n ˜ ), a boundary condition is set which is either metallic (on ∂Ωm , ~n × E ~ ~ ¯ or absorbing (on ∂Ωa , ~n × E = −cµ ~n × (~n × H), where the medium is assumed isotropic, i.e. ε = I3 , ¯ = µI3 and the local light speed c is given by µc2 = 1). µ We assume we dispose of a partition of the domain Ω into a finite number of polyhedral elements. The reader can refer to (Piperno & Fezoui 2003) for a more detailed description of the method in the general case of arbitrary elements. We will focus in the next sections of this paper on orthogonal grids made of orthogonal hexahedra, in particular on partitions which are non-conforming in the sense of Figure 1, where some hexahedra have been regularly divided. ¯i are In the most general case, for each polyhedral element Ti , Vi denotes its volume, and ε¯i and µ respectively the local electric permittivity and magnetic permeability tensors of the medium, which are assumed constant inside the element Ti . We call interface between two neighboring elements their intersection (if it is a surface). One should notice that this interface can be only aT part of a face of one hexahedron in the non-conforming case. For each internal interface aik = Ti Tk , we denote by ~nik the integral over the interface of the unitary normal, oriented from Ti towards Tk . TheSsame definitions are extended to boundary interfaces (in the intersection of the domain boundary ∂Ω m ∂Ωa with an element), the index k corresponding to a fictitious element outside the domain. We denote by t ~n ˜ ik = ~nik /k~nik k = (˜ nikx , n ˜ iky , n ˜ ikz ) the normalized normals and by Vi the set of elements P neighboring Ti (i.e. having an interface in common with Ti ). Let us define the perimeter Pi by Pi = k∈Vi k~nik k.

Figure 1: A typical, orthogonal, locally refined, non-conforming mesh. Inside each element Ti , the numerical unknowns of the method are related to the orthogonal (in the sense of the classical L2 scalar product) projection of the electric and magnetic fields on a chosen local basis of vector fields ϕ ~ ij , 1 ≤ j ≤ di , where di denotes the number of local scalar degrees of freedom. Absolutely no continuity is enforced through element boundaries. Dot-multiplying (1) by a given basis

3

~ X. ~ = curl ~ ψ. ~X ~ × X) ~ ψ ~ − div(ψ ~ yields function ϕ ~ ij , integrating over Ti and using the identity curl  Z Z Z Z ~  ∂E  t ~ ~j.~ ~ ~ ~  ¯ ϕ ~ ij εi = curl ϕ ~ ij .H − (~ ϕij × H).n ˜− ϕij ,  ∂t Ti Ti ∂Ti Ti Z Z Z ~  ∂H t  ~ ϕ ~+ ~ ~n  ¯i =− curl ~ ij .E (~ ϕij × E). ˜. ϕ ~ ij µ  ∂t Ti ∂Ti Ti

(2)

~ i and H ~ i , the L2 -orthogonal projections of the fields The equations above can be rewritten in terms of E Ti ~ ~ E and H on the local field space Span(~ ϕij , 1 ≤ j ≤ di ) which are decomposed the following way: X X ~ i (x, t) = ~ i (x, t) = ∀x in Ti , E Eij (t) ϕ ~ ij (x), H Hij (t) ϕ ~ ij (x). (3) 1≤j≤di

1≤j≤di

This rewriting is straightforward for all volume integrals. For boundary integrals, since no continuity is imposed on the fields, some additional approximations have to be done. We choose here to use a centered flux approximation: ~ i (x) + H ~ k (x) H ~ H(x) ' . 2

~ i (x) + E ~ k (x) E ~ , k ∈ Vi , ∀x ∈ aik , E(x) ' 2

Inside each element, the electric and magnetic fields both depend on di scalar values, denoted by Ei = (Eil )1≤l≤di and Hi = (Hil )1≤l≤di . For the time discretization, we use an explicit leap-frog timescheme, which leads to an energy-conserving, dissipation-free numerical method. In the sequel,the time step ∆t is fixed. Superscripts refer to time stations: unknowns related to the electric field (respectively the magnetic field) are approximated at integer time stations tn = n∆t (respectively half-integer time n+1/2 n stations tn+1/2 = (n + 1/2)∆t) and are denoted by Eij (respectively by Hij ). All definitions for Eni , n +1 / 2 n+1/2 ~ n ~ are similarly extended. Finally, (2) is discretized in time and space according Hi , Ei , and H i to   Z Z n+1 XZ ~ n+1/2 + H ~ n+1/2  H − Eni n+1/2  Ei n+1/2 i  k ~ ~ ~  M = curl ϕ ~ . H − j .~ ϕ − (~ ϕ × ).~n ˜ ik , ij ij ij  i i  ∆t 2 T T a j i i ik k∈Vi  (4)  Z n+3/2 ~ n+1 ~ n+1 + E XZ  E − Hin+1/2  M µ Hi n+1 i k ~ ~ ~  (~ ϕij × ~ ij .Ei + = − curl ϕ ).n ˜ ik .  i  ∆t 2 a T j

i

k∈Vi

ik

~ i and H ~i where the j subscripts denote the jth component of vectors in the left hand side, the fields E are given in (3) in function of scalar degrees of freedom, and Mi and Miµ are local square mass matrices of size di , given by Z (Mi )jl

=

(Miµ )jl

=

ZT i

Ti

t

ϕ ~ ij ε¯i ϕ ~ il , 1 ≤ j, l ≤ di ,

t

¯i ϕ ϕ ~ ij µ ~ il , 1 ≤ j, l ≤ di .

It is clear that the matrices Mi and Miµ are symmetric and positive definite, because the tensors ε¯i and ¯i are symmetric positive definite, and the basis functions ϕ µ ~ ij are assumed to be linearly independent. Metallic and absorbing conditions are dealt with in a weak sense by taking some values for the ~ and H ~ inside the fictitious element Tk (neighboring element Ti ), located beyond the boundary fields E ~ n = −E ~ n and H ~ n+1/2 (x) = H ~ n+1/2 (x) on aik . If aik is face aik considered. If aik is metallic, we use E i k i k an absorbing boundary face, numerical fluxes are derived from a first-order Silver-M¨ uller absorbing condition (see (Piperno & Fezoui 2003) for more details). In that case, fictitious fields are taken as ~ n+1/2 , where i , µi and ci are respectively the ~ n and E ~ n+1 (x) = −ci µi ~n ~ n+1/2 = ci i ~n ˜ ik × H ˜ ik × E H i i k k local electric permeability, magnetic permittivity, and light speed (the medium is assumed isotropic and homogeneous near the absorbing boundary). A stability result has been proved in the general case of arbitrary polyhedral grids (Piperno & Fezoui 2003). More precisely, a discrete (modified near absorbing boundaries) equivalent of the electromagnetic energy is proved to be non-increasing under some CFL-like stability condition on ∆t. Also, 4

this energy is exactly conserved if only metallic conditions are used, which proves the scheme is not dissipative. Concerning convergence, the global framework of DGTD methods yields proofs of convergence properties, in some (not much) simplified case where the physical material considered is piecewise constant and local functional sets are uniformly fixed (Cockburn et al. 2000, Cockburn et al. 2004, Fezoui, Lohrengel, Lanteri & Piperno 2005), which is not exactly the case in this study where the grid is non-conformingly locally refined. Perfectly Matched Layers (PML), introduced by Berenger (Berenger 1996), can be used to deal efficiently with the solution of electromagnetic problems in unbounded regions. This method is based on a split-field formulation of Maxwell’s equations. It is now well known that this formulation is weakly well posed (with possible linear growth with time) (B´ecache & Joly 2002). The 12 PML equations can then be discretized similarly, through dot-multiplication by each local basis function ϕ ~ ij , integration over element Ti and centered numerical fluxes. However, an exponential time scheme which includes PML conductivities is used. The absorbing medium itself has to be bounded. We have numerically observed that using a metallic boundary condition on the PML reconstructed fields can lead in some cases to numerical instabilities. These instabilities never appeared when the first order boundary absorbing condition was used.

A particular DGTD method on orthogonal hexahedra We introduce two sets of basis vector fields very well adapted to Cartesian grids (orthogonal hexahedra). More precisely, we use divergence-free polynomial (P1 or P2 ) basis vector fields. For each hexahedron Ti , let us introduce Gi = t(xGi , yGi , zGi ) its mass center and ∆ = t(∆xi , ∆yi , ∆zi ) its sizes. The DG-P1div method defined by (4)-(5) is based on the P1div space generated by the 9 following L2 -orthogonal basis functions:  ~ i1 = t(1, 0, 0), ϕ ~ i2 = t(y − yGi , 0, 0), ϕ ~ i3 = t(z − zGi , 0, 0),  ϕ t t ϕ ~ i4 = (0, 1, 0), ϕ ~ i5 = (0, x − xGi , 0), ϕ ~ i6 = t(0, z − zGi , 0), (5)  t t ϕ ~ i7 = (0, 0, 1), ϕ ~ i8 = (0, 0, x − xGi ), ϕ ~ i9 = t(0, 0, y − yGi ).

We also introduce two other Discontinuous Galerkin methods based on extended spaces of vector fields. We first introduce the space D which is L2 -orthogonally supplemented from P1div by the 3 additional L2 -orthogonal basis functions:  ~ i10 = t((y − yGi )(z − zGi ), 0, 0),  ϕ ϕ ~ i11 = t(0, (x − xGi )(z − zGi ), 0), (6)  ϕ ~ i12 = t(0, 0, (x − xGi )(y − yGi )). The Discontinuous Galerkin DG-D method defined by (4)-(5)-(6) has not been implemented. It is introduced here in order to prove some divergence preservation properties. We finally introduce the space P2div which is L2 -orthogonally supplemented from D by the 6 additional L2 -orthogonal basis functions:  ~ i13 = t((y − yGi )2 − ∆yi2 /12, 0, 0), ϕ ~ i14 = t((z − zGi )2 − ∆zi2 /12, 0, 0),  ϕ ~ i16 = t(0, (z − zGi )2 − ∆zi2 /12, 0), ϕ ~ i15 = t(0, (x − xGi )2 − ∆x2i /12, 0), ϕ (7)  ~ i18 = t(0, 0, (y − yGi )2 − ∆yi2 /12). ϕ ~ i17 = t(0, 0, (x − xGi )2 − ∆x2i /12), ϕ

The DG-P2div method, defined by (4)-(5)-(6)-(7), will be used in this work only in the context of nonconforming locally refined grids. Since the bases of P1div , D, and P2div are L2 -orthogonal, the schemes will only require the ”inversion” of diagonal local mass matrices. Another advantage is that the orthogonal form of the bases makes the hybridation of schemes (like P1div /P2div ) very easy to implement. One can also notice that the spaces Pndiv (n = 1, 2) require a smaller number of degrees of freedom (respectively 9 and 18 for n = 1 and n = 2) than the full spaces Pn (respectively 12 and 30 for n = 1 and n√= 2). Concerning stability, the numerically observed maximal values of the Courant number ν = c∆t 3/h leading to a stable method on cubes of size h can be compared to the corresponding value for Yee’s scheme, for which ν = 1. We have observed that νP1div = 0.65 and νP2div = 0.4. A general sufficient stability condition can be proved for general non-cubic Cartesian grids, but it would be more restrictive. 5

Dispersion analysis. In order to estimate the accuracy of the DG-P1div scheme, we propose a dispersion analysis on uniform cubic grids in the homogeneous isotropic case and compare it to the Yee’s scheme as a reference scheme on orthogonal grids. A general dispersion analysis for DGTD methods on Cartesian grids for complete sets of polynomial functional sets has been given by Ainsworth (Ainsworth 2004). In the present case, the polynomial functional sets are not complete. Because of the difficulty of formally evaluating the eigenvalues of an 18 × 18 matrix (even with the help of available formal mathematical software!), we will restrict the study to plane waves with a wave vector k parallel to one of the principal directions t (1, 0, 0), t (1, 1, 0) or t (1, 1, 1). Details will be given for the first direction (i.e. k =t (1, 0, 0)). In that case, the problem reduces to a one-dimensional system and the unknowns are the components Hy and Ez with two degrees of freedom each. On the element [j∆x, (j + 1)∆x], let us set Vj =t (Hj,1 , Hj,2 , Ej,1 , Ej,2 ), the indices y and z being omitted for H and E respectively. √ As a plane wave, Vj takes the form Vj = V exp(iω(n + s)∆t + ijk∆x) where V is the amplitude, i = −1 ,k = |k|, s = 0 for the electric components and s = 12 for magnetic ones. By injecting the plane wave Vj into the DG-P1div  2i scheme, we obtain that λ = c∆t sin ω∆t is one of the eigenvalues of a 4 × 4 matrix (see (Canouet 2003) 2 for details). Using Taylor expansions, one finds that two eigenvalues lead to propagative waves while √ the two others correspond to oscillating spurious modes. In function of ν = c∆t 3/∆x and the exact pulsation ωex = kc, one finds the following dispersion relations for each principal direction: 2

t

2 k k (1, 0, 0) : ω 2 /ωex = 1 + k 2 ∆x2 ( ν36 + 2 t 2 2 k k (1, 1, 0) : ω /ωex = 1 + k 2 ∆x2 ( ν36 − 2 t 2 k k (1, 1, 1) : ω 2 /ωex = 1 + k 2 ∆x2 ( ν36 −

1 4 4 24 ) + O(k ∆x ), 1 4 4 168 ) + O(k ∆x ), 1 4 4 36 ) + O(k ∆x ).

We recall the corresponding dispersive relations of Yee’s scheme: 2

t

2 k k (1, 0, 0) : ω 2 /ωex = 1 + k 2 ∆x2 ( ν36 − 2 t 2 2 k k (1, 1, 0) : ω /ωex = 1 + k 2 ∆x2 ( ν36 − 2 t 2 k k (1, 1, 1) : ω 2 /ωex = 1 + k 2 ∆x2 ( ν36 −

1 12 ) 1 24 ) 1 36 )

+ O(k 4 ∆x4 ), + O(k 4 ∆x4 ), + O(k 4 ∆x4 ).

We show in Figure 2 the second order term of the dispersion error as a function of the Courant number 1

0.8

Yee-100

ERR

DG-100

0.6

0.4 Yee-110 DG-111

0.2

Yee-111 DG-110

0 0

0.2

0.4

NU

0.6

0.8

1

Figure 2: DG-P1div scheme and Yee’s scheme: 2nd order dispersion term for principal k. ν for both the DG-P1div scheme (ν ∈ [0, 0.65]) and Yee’s scheme (ν ∈ [0, 1]). The results are normalized 2 using the maximum error of Yee’s scheme, i.e. we plot err = 12 |ω 2 /ωex − 1|/(k 2 ∆x2 ) in function of ν. Although this study is incomplete because it is only based on observations concerning three principal directions of propagation, we can make a few remarks. First, near the limit of the stability of each scheme, the dispersion errors of both schemes are really close to each other. This means that the 6

performance of both schemes will be comparable when the maximal admissible ν is used. Second, as ν gets smaller, the maximal dispersion increases for Yee’s scheme while it is reduced for the DG-P 1div scheme (when ν → 0, its minimal dispersion error is really small and its maximal dispersion is close to the minimal dispersion error of Yee’s scheme). This means that using small time steps will not damage the precision of the DG-P1div scheme. Thus there is no mandatory need to use local time stepping when Maxwell’s equations are solved on locally refined grids, which is not the case for Yee’s scheme. Divergence preservation. In the particular case of Cartesian grids, it can be shown that the DGTD ~ = ρ, method proposed preserves in some weak sense the Gauss divergence laws given by div( ε¯E) ∂ρ ~ ~ ¯H) = 0, where the density of electric charge ρ verify ∂t = div(j). It is clear that Gauss laws div(µ are automatically satisfied by transient solutions of Maxwell’s equations if they are satisfied at initial time. However, since Maxwell’s equations are only solved in an approximate way, it is very interesting to investigate whether numerical approximations of Gauss laws are verified by the numerical approximate fields. We will prove in this section that the DG-D method preserves the Gauss laws in a weak sense and that the DG-P1div method can be viewed as a restriction, in terms of degrees of freedom, of the DG-D method. We define the classical Q1 space of scalar functions which are continuous, with compact support included in Ωh , and polynomial of degree at most one in each space coordinate. We first prove the following theorem: Theorem 1 (Weak preservation of Gauss laws) We consider the Discontinuous Galerkin DG-D method defined by (4)-(5)-(6) on a conforming orthogonal hexahedral grid Ω h . Then, the approximate ~ n, H ~ n+1/2 )i,n verify electromagnetic fields (E i

i

Z  Z ~ n+1 − E ~ n ) = ∆t  ¯(E ψ div~jn+1/2 , ψ div ε   Ω Ω h h ∀ψ ∈ Q1 , Z   n +1/2 n -1/2 ~ ~  ¯ (H ψ div µ −H ) = 0.

(8)

Ωh

Proof. It is sufficient to prove (8) for any Q1 basis function ψ naturally linked to a given vertex S inside the conforming orthogonal hexahedral grid. Let us define the support of ψ as C = ~ ∇ψ ~ lies in D, that curl ~ = ~0 and {Ti , S is a vertex of Ti }. It is very important to notice that ∇ψ ~ that since ψ is continuous, the tangential jump of ∇ψ through faces of neighboring hexahedra is zero, ~ ik × ~n i.e. [∇ψ] ˜ ik = ~0 for any rectangular face between adjacent hexahedra. We recall we denote by T aik = Ti Tk the interface between neighboring hexahedra. We have successively: Z Z Z n+1/2 n−1/2 n+1/2 n−1/2 ~ ~ ~ ~ ~ n+1/2 − H ~ n−1/2 ) ¯ ¯ ¯ (H ψ divµ(H −H )= ψ divµ(H −H )=− ∇ψ.µ Ωh C C !! XZ X Z ~n +E ~n E n i k ~ ~ × ~n ˜ ik curl(∇ψ = ∆t ∇ψ|Ti . |Ti ).Ei − 2 Ti k∈Vi aik Ti ⊂C !! X XZ ~n +E ~n E i k × ~n ˜ ik = −∆t ∇ψ|Ti . 2 Ti ⊂C k∈Vi aik !! Z   X ~n+E ~n E i k ~ |T − ∇ψ ~ |T × ~n = ∆t = 0. ∇ψ ˜ ik . i k 2 aik aik ⊂C/(∂C)

R R ~ n+1 − E ~ n ) = ∆t ψ div~jn+1/2 , which concludes the proof.  Similarly, we obtain Ωhψ divε¯(E Ωh One can notice that the DG-P1div scheme does not verify (8). Nevertheless, it is easy to check that the methods DG-P1div and DG-D are partly identical. More precisely, for all elements Ti and all faces aik , we have: Z Z Z ~ ~ ∀ 1 ≤ j ≤ 9, ∀ 10 ≤ l ≤ 12, curl(~ ϕij ).~ ϕil = 0, (~ ϕij × ϕ ~ il ).n ˜ ik = 0, (~ ϕij × ϕ ~ kl ).~n ˜ ik = 0. Vi

aik

7

aik

15

15

Hy

10

10

5

5

0

0

-5

-5

-10

-10

-15

0

20

40

60 80 100 Time (ns)

120

140

-15

160

Hy

0

20

40

60 80 100 Time (ns)

120

140

160

Figure 3: P1div method (left) and P1div /P2div method (right): Hy (t) inside the refined grid (ratio 8). As a consequence, one can verify that the first 9 components of the fields in each element computed by the two schemes are identical. Then the DG-P1div scheme can be seen as a computationally more efficient restriction of the DG-D scheme, which preserves divergence laws in a weak sense.

Elementary tests and verification of numerical properties Using the P1div approximation with a conforming refinement is very efficient. For non-conforming grids however, spurious reflections could be produced by non-conforming interfaces. We consider a radiating dipole located at the center of a [0, 1]3 domain in the free space. The domain is discretized using a 21×21×21 cubic grid, complemented by a eight cell PML layer. It is locally refined in a non-conforming way in the zone [0.24, 0.43]3 with a refinement ratio equal to 8. For this experiment, the source current ~j R is parallel to the z axis, located inside a single hexahedron and computed using Ti ~j.~ ϕij = δj7 ∂t ρ, where 2

ρ(t) = 10−10 exp (−(t/T − 3) ) (with T = 2ns). When the DG-P1div method is used, some spurious reflections appear and some energy is also captured inside the refined region. We have checked that these reflections and the captured energy are small compared to the original signal. Anyway, we show in Figure 3 (left) the component Hy computed using the P1div method, in function of the time at a point located inside the refined grid. Spurious undamped oscillations are present after the pulse is gone through the refined zone. One can conclude that the P1div method is unfortunately unable to handle accurately non-conforming grids. A first strategy consists in increasing the order of the approximation using the P 2div approximation for example. This approach would require a larger number of degrees of freedom, and additional CPU time and memory costs. We propose here to locally use the P2div approach only in elements of the coarse grid which neighbor the refined grid. Everywhere else, the P1div approximation is still used, as described in Figure 4. We name the resulting scheme the DG-P1div /P2div method. The stability of such hybrid scheme is ensured under some CFL-like condition. It was numerically observed that the hybrid scheme has the same stability limit (ν ≤ 0.65 based on refined cubes) as the DG-P1div method, as long as the grid is actually refined. This is not a surprise, since the DG-P2div method, which has a reduced stability domain (somewhere around ν ≤ 0.4 on regular cubes) is only used on some elements of the coarse grid (which are at least twice larger than elements of the refined grid). The solution provided by the DG-P1div /P2div method on the same test-case is proposed on the right part of Figure 3. The spurious oscillations have disappeared. This is also illustrated on Hx contours on the z = 0.3 cut plane for several times, in Figures 5 and 6, respectively for the DG-P1div method and the DG-P1div /P2div method. The rectangle shows the limit of the refined zone. The effects of the numerical dispersion should also be evaluated on locally refined grids. In a one meter cubic resonant cavity, we use the DG-P1div /P2div method to compute the time evolution of the (1, 1, 1)-mode (standing wave at 0.26 GHz) using a non-conforming locally refined cubic grid with ratio

8



1

Pdiv













                                                                                                                                                                                                                                                                                                                                                                                                                                                                          

                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               

                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         

2

1

Pdiv

Pdiv

Figure 4: P1div /P2div : two-dimensional non-conforming example.

Figure 5: GD-P1div : Hx component in cut plane z = 0.3 at times 4ns, 4.7 ns, 5.3ns, 6ns, 6.5ns, 7ns.

9

Figure 6: GD-P1div /P2div : Hx component in cut plane z = 0.3 at times 4ns, 4.7 ns, 5.3ns, 6ns, 6.5ns, 7ns. 10 (i.e. ratio 1000 for cube volumes). The number of points per wavelength λ/∆x is equal to 15 in the coarse grid and 150 in the refined zone. The degrees of freedom are initialized by the projection of the exact solution on the local basis functions. We compare in Figure 7 the approximate and exact 0.4 Discontinuous Galerkin Exact

0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 45

46

47

48

49

50

Figure 7: Ez (t, x0 ) - zoom on the last 5 periods - refinement ratio 10. components Ez in function of the time at a given point in the coarse grid (the 5 last periods of 45 are shown). The curves are very close. This shows the scheme does not introduce too much dispersion in the coarse grid where ν is far from optimal, although the refinement ratio is high and the simulation is quite long. We show in Figure 8 the component Ez along a diagonal cut plane for the exact and the approximate solutions (after projection on the coarse grid). One can see that they compare very well.

10

Discontinuous Galerkin

Exact solution

Figure 8: Approximate (left) and exact (right) Ez fields (plane cut x + y + z = 3/2). In order to evaluate the influence of the non-conforming refinement ratio, numerical experiments for the same cavity, using the DG-P1div /P2div method have been made with several refinement ratios with λ/∆x = 15 in the coarse grid. It appears that the L2 norm of the difference between the approximate and exact solutions does not depend on the refinement ratio (see Figure 9). 0.16

Refinement 1:2 Refinement 1:4 Refinement 1:6 Refinement 1:8

0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0

5

10

15 Time

20

25

30

Figure 9: L2 error on cavity mode simulations for several refinement ratios.

Numerical simulations in industrial configurations Patch antenna. We first study the resonant frequency of the printed antenna described in Figure 10. The a×a square patch and the l×b microstrip line are printed above a L×L metallic ground plane. Both planes are separated by a fictitious dielectric layer (thickness h, relative electric permittivity ε r = 1). A theoretical approximate for the resonant frequency f of the antenna is f =5.6 GHz. The electric source is a vertical dipole located at the end of the line (at h/2 from both planes). The transient amplitude of the dipole is Gaussian (frequency 5.6 GHz and bandwidth 2 GHz). The simulation is stopped after 23 ns, when the overall energy is close to zero. A Fourier analysis of fields between the patch and the ground plane yields the approximate resonant frequency, which depends on the grid used. In all cases, a surrounding 6-cell thick PML has been used and the corresponding PML elements are included in the total numbers of elements mentioned. Simulations on three different grids are performed. In 11

Ground plane a

Patch a = 26.8 mm b = 6.69 mm b

L

h = 2.23 mm l = 42.4 mm

Line

L = 80.4 mm

l

Dipole Figure 10: Square patch fed by a microstrip line. all cases, the dimensions of the computational domain are 147.4mm × 147.4mm × 67mm. The locally refined grids considered are located in a 40.14mm × 89.2mm × 8.92mm zone which encloses the whole line, the square patch and a part of the ground plane. The main characteristics of these grids and the approximate frequencies obtained are given in Table 1. It clearly appears that geometrical details can Table 1: Characteristics of grids used for the patch antenna. Grid number Grid type Refinement ratio λ/∆x (coarse grid) λ/∆x (refined grid) # elements equivalent # elements ∆xmin ∆xmax Elements between planes Approximate frequency

#1 uniform 24 24 130,680 130,680 h h 1 5.33 GHz

#2 locally refined 4 12 48 52,301 1,045,440 h/2 2h 2 5.46 GHz

#3 locally refined 6 12 72 107,021 3,528,360 h/3 2h 3 5.54 GHz

be taken into account more accurately, and then yield a more accurate resonant frequency (5.54 GHz is to be compared with the exact frequency 5.56 GHz), when the grid is refined. The local non-conforming refinement makes it possible to have a finer grid wherever it is necessary, with a reduced total number of elements required. We have represented in Figure 11 contours of Ez at t=10 ns on metallic elements of the non-conforming locally refined grid #3 (the distance between planes has been magnified).

Printed dipole with a slot. Finally, we model in this last example the printed dipole with a slot presented in (Brachat, Dedeban & Ratajczak 1997). The central working frequency of this structure is 1.15 GHz. It is made of printed elements on both sides of a thin dielectric λ × λ/2 rectangular plate, with ε r = 2.2 and thickness h = λ/160. On the upper side of the dielectric, one has printed a ”T” with a slot of width f = λ/500. On the lower side, a feeding line is printed. The feeding source is a dipole located at the end of the line. It is radiating with a central frequency of 1 GHz and a half bandwidth of 0.5 GHz. We perform a simulation of this model using a non-conforming locally refined grid. In order to be as accurate as possible (fine grid, absorbing layers as far as possible), we have chosen to take in the fine grid ∆x = λ/500 (which yields one element inside the slot), ∆y = λ/250, and ∆z = λ/333 (which yields two 12

Figure 11: Ez at t = 10 ns on metallic faces of the locally refined grid 3 with ratio 6. elements in the dielectric thickness). In the coarse grid, we take ∆x = ∆y = ∆z = λ/42, which yields different refinement ratios in the three directions (respectively 12, 6, and 8). The absorbing boundary is set at λ/4 (x and y directions) or λ/3 (z direction). The global grid has 571,000 elements. The equivalent uniform cubic grid would be made of 93 million elements (with ∆x = ∆y = ∆z = λ/500). ~ on metallic faces (”T” and line) at 900 MHz. They are presented We compute electric currents ~n × H in Figure 12. For the same structure, similar currents are shown in Figure 13. They were obtained using the SR3D software of France Telecom R&D (Brachat et al. 1997, Ratajczak, Brachat & Guiraud 1994) (based on integral equations (Bendali 1984, N´ed´elec 1980)) with a 20,000 triangular, unstructured, surfacic mesh of metallic elements. The results obtained compare very well. In both cases, currents are almost constant along the line and maximal currents are obtained near the end of the slot. For this simulation, the optimized, parallel software SR3D requires a CPU time of one hour (on 6 hppe-8007 750 MHz processors). Then a broadband study would require at least 15 similar computations at 15 different frequencies (between 800 MHz and 1.5 GHz, with a step of 50 MHz). The same study could be done in one time-domain simulation. The present time-domain simulation based on a sequential implementation required a CPU time of 70 hours (on a 2.4 GHz Pentium 4 PC), which would be reduced after parallelization. This makes the DGTD volumic approach competitive. We can also notice that the cost of the DGTD approach does not depend on material discontinuities, which is not the case for integral equation codes. This makes the volumic DGTD approach even more interesting for highly heterogeneous structures. In order to put the emphasis on the advantages of locally refined grids, we propose a computation in a similar configuration with a transparent dielectric material (i.e. εr = 1). In that case, the refinement zone can be concentrated along the metallic elements of the feeding line. We use ∆x = ∆y = ∆z = λ/500 in the fine grid and ∆x = ∆y = ∆z = λ/45 in the coarse grid (refinement ratio of 11 in all three directions). The global grid is made of 655,000 elements and the equivalent cubic uniform grid would be made of 86 million elements. The currents are presented in Figure 14 for the frequency 900 MHz. The contours have changed but seem coherent. One can notice that the slot is well simulated, although the non-conforming interface has been set closer. At the same time, the currents obtained on metallic

13

Figure 12: Currents at 900 MHz (printed dipole with dielectric material).

Figure 13: SR3D software: mesh and currents at 900 MHz (printed dipole with dielectric material).

14

Figure 14: P1div /P2div DGTD method: currents at 900 MHz (printed dipole with no material). faces which are inside the coarse grid seem also quite accurate.

Conclusion We have described a fully explicit Discontinuous Galerkin Time Domain (DGTD) method which has been applied to solve Maxwell’s equations on both conforming and non-conforming grids. The method conserves a discrete energy and is stable under a CFL condition. It also preserves divergence relations in a weak sense. Finally, numerical dispersion errors do not increase significantly when using small time steps. The numerical experiments we carried out revealed the high quality of the results together with the flexibility of the method, allowing us to mix different sets of local basis functions, to use PML or upwinding techniques at absorbing boundaries, and finally to use highly refined grids without a noticeable loss of accuracy compared to uniformly refined grids.

Acknowledgments The authors thank C. Dedeban for his priceless help on many steps along this work and P. Ratajczak for his assistance on the SR3D software of France Telecom R&D.

References Ainsworth, M. (2004), ‘Dispersive and dissipative behaviour of high order discontinuous Galerkin finite element methods’, J. Comput. Phys. 198, 106–130. B´ecache, E. & Joly, P. (2002), ‘On the analysis of B´erenger’s perfectly matched layers for Maxwell’s equations’, RAIRO Mod´el. Math. Anal. Num´er. 36(1), 87–119.

15

Bendali, A. (1984), ‘Numerical analysis of the exterior boundary value problem for the time harmonic Maxwell equations by a boundary finite element method’, ”Math. Comp.” 43, 29–46. Berenger, J.-P. (1996), ‘Three-Dimensionnal Perfectly Matched Layer for the Absorption of Electromagnetic Waves’, J. Comput. Phys. 127, 363–379. Brachat, P., Dedeban, C. & Ratajczak, P. (1997), Modeling of arbitrary radiating elements using the CAD tool SR3D, in ‘IEEE Antennas and Propagation Symposium’, Montr´eal, Canada, pp. 1802– 1805. Canouet, N. (2003), M´ethode de Galerkin Discontinu pour la r´esolution du syst`eme de Maxwell sur des ´ maillages localement raffin´es non-conformes, math´ematiques appliqu´ees, Ecole Nationale des Ponts et Chauss´ees. Cioni, J.-P., Fezoui, L., Anne, L. & Poupaud, F. (1997), A parallel FVTD Maxwell solver using 3D unstructured meshes, in ‘13th annual review of progress in applied computational electromagnetics’, Monterey, California, pp. 359–365. Cockburn, B., Karniadakis, G. E. & Shu, C.-W., eds (2000), Discontinuous Galerkin methods. Theory, computation and applications., Vol. 11 of Lecture Notes in Computational Science and Engineering, Springer-Verlag, Berlin. Cockburn, B., Li, F. & Shu, C.-W. (2004), ‘Locally divergence-free discontinuous galerkin methods for the Maxwell equations’, J. Comput. Phys. 194(2), 588–610. ´ ements finis d’arˆete et condensation de masse pour les ´equations de Elmkies, A. & Joly, P. (1997), ‘El´ Maxwell: le cas de dimension 3.’, C. R. Acad. Sci. Paris S´er. I Math. t. 325(11), 1217–1222. Fezoui, L., Lohrengel, S., Lanteri, S. & Piperno, S. (2005), ‘Convergence and stability of a Discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes’, submitted for publication to RAIRO Mod´el. Math. Anal. Num´er. Hesthaven, J. & Teng, C. (2000), ‘Stable spectral methods on tetrahedral elements’, SIAM J. Sci. Comput. 21(6), 2352–2380. Hesthaven, J. & Warburton, T. (2002), ‘Nodal high-order methods on unstructured grids. I: Timedomain solution of Maxwell’s equations.’, J. Comput. Phys. 181(1), 186–221. Hesthaven, J. & Warburton, T. (2004), ‘High-order nodal discontinuous Galerkin methods for the Maxwell eigenvalue problem’, Philos. Trans. Roy. Soc. London Ser. A 362, 493–524. Holland, R., Cable, V. P. & Wilson, L. C. (1991), ‘Finite-volume time-domain (FVTD) techniques for EM scattering’, IEEE Trans. Electromagn. Compat. 33(4), 281–294. Hyman, J. M. & Shashkov, M. (1999), ‘Mimetic discretizations for Maxwell’s equations’, J. Comput. Phys. 151, 881–909. Joly, P. & Poirier, C. (2000), A new second order 3D edge element on tetrahedra for time dependent Maxwell’s equations, in A. Bermudez, D. Gomez, C. Hazard, P. Joly & J.-E. Roberts, eds, ‘Fifth International Conference on Mathematical and Numerical Aspects of Wave Propagation’, SIAM, Santiago de Compostella, Spain, pp. 842–847. Kopriva, D. A., Woodruff, S. L. & Hussaini, M. Y. (2000), Discontinuous spectral element approximation of Maxwell’s equations, in B. Cockburn, G. E. Karniadakis & C.-W. Shu, eds, ‘Discontinuous Galerkin methods. Theory, computation and applications.’, Vol. 11 of Lecture Notes in Computational Science and Engineering, Springer-Verlag, Berlin, pp. 355–362. Lee, J.-F. & Sacks, Z. (1995), ‘Whitney elements time domain (WETD) methods’, IEEE Trans. Magn. 31(3), 1325–1329.

16

Lohrengel, S. & Remaki, M. (2002), A FV Scheme for Maxwell’s equations: Convergence Analysis on unstructured meshes, in R. Herbin & D. Krner, eds, ‘Finite Volumes for Complex Applications III’, Hermes Penton Science, London, Porquerolles, France, pp. 219–226. N´ed´elec, J.-C. (1980), ‘Mixed finite elements in R3 ’, Numer. Math. 35, 315–441. Piperno, S. & Fezoui, L. (2003), A Discontinous Galerkin FVTD method for 3D Maxwell equations, Technical Report RR-4733, INRIA. Piperno, S., Remaki, M. & Fezoui, L. (2002), ‘A non-diffusive finite volume scheme for the 3D Maxwell equations on unstructured meshes’, SIAM J. Numer. Anal. 39(6), 2089–2108. Ratajczak, P., Brachat, P. & Guiraud, J.-L. (1994), ‘Rigorous analysis of 3D structures incorporing dielectrics’, IEEE Trans. Antennas and Propagation 42(8). Remaki, M. (2000), ‘A new finite volume scheme for solving Maxwell’s system’, COMPEL 19(3), 913– 931. Remaki, M. & Fezoui, L. (1998), Une mthode de Galerkin Discontinu pour la rsolution des quations de Maxwell en milieu htrogne, Technical Report RR-3501, INRIA. Warburton, T. (2000), Application of the discontinuous Galerkin method to Maxwell’s equations using unstructured polymorphic hp-finite elements., in B. Cockburn, G. E. Karniadakis & C.-W. Shu, eds, ‘Discontinuous Galerkin methods. Theory, computation and applications.’, Vol. 11 of Lecture Notes in Computational Science and Engineering, Springer-Verlag, Berlin, pp. 451–458. Yee, K. S. (1966), ‘Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media’, IEEE Trans. Antennas and Propagation AP-16, 302–307. Yioultsis, T. V., Kantartzis, N. V., Antonopoulos, C. S. & Tsiboukis, T. D. (1998), ‘A fully explicit Whitney element time domain scheme with higher order vector finite elements for three-dimensional high frequency problems’, IEEE Trans. Magn. 34(5), 3288–3291.

17