A  novel  implicit  nodal  Space-Time  Discontinuous
Galerkin (STDG) method is proposed in this paper. An eigenvalue
analysis is performed and compared with that for a DG scheme
solved  with  a  4th-Order  Runge-Kutta  time  integrator.  We  show
that  STDG  offers  a  significant  improvement  of  dissipative  and
dispersive  properties  and  allows  larger  time  steps,  regardless  of
the  spatial  hp-refinement.  A  domain-decomposition  technique  is
used to introduce an explicit formulation of the method in order
to  render  it  computationally  efficient.This  work  is  partially  funded  by  the  National  Projects  TEC2010-20841-
C04-04,  TEC2013-48414-C3-01,  CSD2008-00068,  P09-TIC-5327,  P12-TIC-
1442, and from the GENIL excellence network

Alvarez Gonzalez, Jesus

Díaz Angulo, Luis Manuel

Fernández Pantoja, Mario Alberto

González García, Salvador

English

Fernandez Pantoja, Mario

Gonzalez Garcia, Salvador

Diaz Angulo, Luis

Institutional Repository of the University of Granada

An Explicit Nodal Space-Time Discontinuous Galerkin Method for Maxwell’s Equations

Repositorio Institucional Universidad de Granada

1An Explicit Nodal Space-Time DiscontinuousGalerkin Method for Maxwell’s EquationsL. D. Angulo, J. Alvarez, Member, IEEE, M. F. Pantoja, Senior Member, IEEE, and S. G. Garcia, Member, IEEEAbstract—A novel implicit nodal Space-Time DiscontinuousGalerkin (STDG) method is proposed in this paper. An eigenvalueanalysis is performed and compared with that for a DG schemesolved with a 4th-Order Runge-Kutta time integrator. We showthat STDG offers a significant improvement of dissipative anddispersive properties and allows larger time steps, regardless ofthe spatial hp-refinement. A domain-decomposition technique isused to introduce an explicit formulation of the method in orderto render it computationally efficient.Index Terms—Discontinuous Galerkin Methods, Computa-tional Electrodynamics, Space-Time Discontinuous GalerkinI. INTRODUCTIONA common approach for Discontinuous Galerkin (DG) time-domain (TD) methods, is to treat the time and space variablesseparately [1], often using an explicit time-integrating schemessuch as the 4th-Order Low-Storage Runge-Kutta LSERK4,which prevents the full exploitation of the higher-order spatialconvergence. Although works to cope with this limitation[2] exist, a noteworthy alternative is to use Space-Time DG(STDG) methods, already used in other fields of Physics[3], [4], and in Electromagnetics [5], [6]. A major drawbackof STDG resides in its implicit nature, though semi-explicitapproaches also exist [5], [6] (tent-pitching technique).In this work, we present a novel STDG formulation com-bined with a causal domain-decomposition technique [7] torender it explicit (E-STDG). This letter is organized as follows:we first formulate a nodal [1] STDG scheme, with a newspurious-free upwind-in-space flux, combined with a centered-in-time flux. Next, we study the properties of the resultingimplicit STDG scheme with an eigenvalue analysis, comparingwith a DG-LSERK4 one. Next, we describe the explicit causalformulation E-STDG, and validate it with a simple resonantproblem. We finally conclude that the use a nodal approach,together with the domain-decomposition technique providesan affordable solution to the problem.II. IMPLICIT FORMULATIONLet us start by defining a 2D space-time region tessellatedwith k = 1, . . . ,Kx elements. With ζ ∈ V defined as a space-time coordinate within each element, and with n̂ = (nt, nx)This work is partially funded by the National Projects TEC2010-20841-C04-04, TEC2013-48414-C3-01, CSD2008-00068, P09-TIC-5327, P12-TIC-1442, and from the GENIL excellence network.The 1st, 3rd, 4th authors are with the Dept. of Electromagnetism, Univ.Granada, 18071 Granada, Spain, (Corresponding: salva@ugr.es). The 2nd iswith Airbus Defense and Space, 28906 Getafe, Spain.The authors thanks Stephane Lanteri, Sascha Schnepp and Martin Lilienthalfor their comments, discussions and support.normal vectors pointing outwards from its boundary. Theweak form of 1D Maxwell’s curl equations is found bymultiplying by weighting functions αij(ζ), integrating overV , and enforcing the residual to vanish. For instance, the free-space Ampère’s law (taking for simplicity the permittivity andpermeability both equal to 1) becomes∫V(∂tE(ζ) + ∂xH(ζ))αij(ζ)dζ = 0 (1)with E and H being the electric and magnetic fields. Integrat-ing by parts in (1) and replacing the boundary flux-integral bya numerical flux, as usual in DGTD, we find∫V∂xHαijdζ =∫Tdt[∫∂Xn̂ ·H∗αijdx−∫XH∂xαijdx](2)with H∗ (and similarly E∗) being the usual numerical fluxes[1]. Next, following the Galerkin procedure, we expand thefields in series using the weighting functions also as basisfunctions, and assuming that their space-time dependence canbe separated in polynomials of orders P t and P x for thetemporal and spatial parts, respectivelyHh(ζ) =P t∑i=0Px∑j=0αij(ζ)Hij =P t∑i=0Px∑j=0αti(t)αxj (x)Hij (3)and similarly for E. Substituting (3) into (2), the spatial-stiffness term becomes,∫VHh∂xαmndζ =P t∑i=0Px∑j=0Hij∫Tαtiαtmdt∫Xαxj ∂xαxndx (4)And the spatial-flux term,∫T(H∗(t, xR)αmn(t, xR)−H∗(t, xL)αmn(t, xL)) dt =P t∑i=0Px∑j=0H∗ij (αxn(xR)− αxn(xL))∫Tαtiαtmdt (5)Flux evaluations are highly simplified by using a nodalapproach based on Lagrange interpolating polynomials [1],since the flux in a space-time node requires only the degreesof freedom (DOFs) at the nodes occupying the same positionin the neighboring elements. To determine the space-time flux,let us define[[Hij ]]t =Hij −Ht,+ij2[[Hij ]]x =Hij −Hx,+ij2(6)and similarly for the E-field. There, the upper subscript x,+refers to the neighboring node along the spatial boundary (see2Fig. 1: Notation: a space-time element (in gray) with P t = P x = 3.Fig. 1). First, we find the centered version of the numericalfluxes, as the average with the contiguous-space/earlier-timeborder valuesHij = δi0[[H0j ]]t − δ0j [[Hi0]]x + δjPx [[HiPx ]]x (7)and, secondly, we complete them with the lacking terms toobtain upwind fluxes in space, keeping centered fluxes in timeH•ij = Hij − δ0j [[Ei0]]xZ−1 + δjPx [[EiPx ]]xZ−1 (8)where δij refers to the usual Kronecker-delta, and Z, Y themedium impedance and admittance, respectively (both unityin our case).A compact matrix formulation can be written by arrangingthe E and H field coefficients (DOFs) in an ordered column-vector, e.g. E = [E0,0, . . . , E0,Px . . . EP t,0, . . . , EP t,Px ]T .Eq. (4), with the usual tensor product ⊗, becomes∫VHh(ζ)∂xα(ζ)dζ =Mt ⊗ CxH (9)with M and C being the mass and stiffness matrices, respec-tively. The spatial-flux term (5) can be expressed as∫T(Hh,∗(t, xR)α(t, xR)−Hh,∗(t, xL)α(t, xL))dt=Mt ⊗ (RxL −RxR)H∗(10)with RtL = ePt0 eP t,T0 and RtR = ePtP teP t,TP t where eNi is aN+1 long zero vector with 1 in entry i. Operating similarly forthe rest of terms and the Ampere’s equation, we can expressthe scheme locally asεCt ⊗MxE +Mt ⊗ CxH=Mt ⊗ (RxL −RxR)H∗ + (RtL −RtR)⊗MxE∗(11)Let us now assume that our computational domain is dividedby Kx space-time elements and let us defineF tR,L = IKx⊗I2⊗RtR,L⊗Mx , FxR,L = IKx⊗V⊗Mt⊗RxR,LSt = IKx⊗I2⊗Ct⊗Mx , Sx = IKx⊗V⊗Mt⊗Cx (12)with IN =∑N−1n=0 eN−1n eN−1,Tn and V = e10e1,T1 + e11e1,T0 .We can write the following compact expression(St + Sx)q = (F tL −F tR + FxL −FxR)q∗ (13)where q represents all the field coefficients in a given space-time element. The upwind-in-space centered-in-time numericalflux (8) can be expressed asq∗upwind = −EtRq′ + (EtL + ExL − ExR)q︸ ︷︷ ︸q∗centered+(ExL − ExR)Ṽq (14)with q′ being the state vector of the previous space-timeelement. The operators E perform the operations needed toassemble the unknowns associated with the fluxes in the globalsystem of equations. The superscripts t and x and subscriptsL and R indicate the boundary at which they are operating.The operator Ṽ = IKx⊗V⊗IP t⊗IPx indicates the operationon the dual field. Note that FxEtq′ = F tExq = 0 becauseE is defined to extract only the unknowns needed by the fluxacting on the boundary indicated by its superscript. Replacingthe centered part of (14) into (13)(St+ Sx−FxLExL + FxRExR + F tLEtL)−1(−F tLEtR)︸ ︷︷ ︸Hq′= q (15)and a similar expression can be formulated for upwind flux.Eq. (15) enables us to find the current state, q, from theprevious state, q′, in a marching-in-time way. A major lim-itation of this scheme is its spatial implicitness, requiringmatrix inversion (or solving a linear system at each time step)of operators scaling as O(K2x) (only tentatively efficient forproblems where the space and time scales are very dissimilar).III. AN EXPLICIT SCHEMESemi-explicit formulations of STDG exist [5], [6], and inthis section we present a new explicit alternative simpler thanthose, in terms of the needed mesh, and using a reducednumber of DOFs. A domain-decomposition technique [7] willbe used to formulate an explicit variant of the scheme in(15), at the cost of adding a CFL-like causality condition.Rather than solving the whole domain Ω, we divide theproblem into smaller regions Ωd such that⋃Ndd Ωd = Ω.Then, we enlarge each region Ωd to include the elements thatare causally connected with it for a given ht, we will denotethis enlarged region with Ω′d. An element is considered to becausally connected with another one if any of its fields at t canpropagate to Ωd at t+ ht. i. e. if it contains a point within adistance d ≤ cht to any point in Ωd, with c being the numericalspeed of light. The border of Ω′d ends in a zero-flux boundarycondition that decouples it from the rest of the computationaldomain. Finally, the problem is solved in each Ω′d using theimplicit method described above but passing only the valuesin Ωd to the next step. The values in Ω′d\Ωd are discardedas they are corrupted by the artificial boundary condition.This approach reduces the computational cost from O(K2x) to∑Ndd O(K2x,d) with Kx.d < Kx. The region sizes Ωd can bechosen as a trade-off taking into account the computationalcost in determining the initial Hd, the cost to evolve thescheme, and the size of the time step.IV. NUMERICAL ANALYSISA PEC-terminated spatial domain has been discretized withKx = 4 elements of size hx = 0.25 and order P x = 4using the spatial upwind flux (8), as a simple proof-of-concepttest-case. The eigenvalues of H have been numerically foundto investigate dispersion and dissipation independently (ratherthan using error norms such as in [3]–[6]). We have conductedstudies for different ht, which have not exhibited eigenvalues3with positive real parts, thus implying that the scheme isunconditionally stable for any ht.Fig. 2 shows the dispersion and dissipation properties of thescheme for the first resonant modes to the closest analyticalmode ka = π with ωa = π for different STDG schemes anda classical DG evolved using LSERK4 and the same upwindspatial fluxes. The accuracy of the scheme presents high-orderconvergence with ht, following the relationship (h2Pt+1t ) forthe dissipation, and (h2Pt+2t ) for the dispersion, in agreementwith [8] for the spatial DG semi-discrete scheme. Thereforewe can conclude that the spatial and temporal convergences ofthe scheme coincide. Note that, although the LSERK4 physicaleigenvalue can be computed and represented, it is unstable forapprox. ht ≥ 0.05 because of the presence of other eigenvalueslying outside its stability region, thus limiting with ht,max ∝(P x)−2 minhkx. We also find that we have higher convergencefor P t ≥ 2, which is a significant improvement over LSERK4,especially when combined with a higher P x.10−2 10−1 10010−1010−810−610−410−2100ht−real(log(H))STDG−Pt 1STDG−Pt 2STDG−Pt 3STDG−Pt 4LSERK4O(h )tO(h )tO(h )tO(h )t357910−2 10−1 10010−1010−810−610−410−2100ht1−imag(log(H))STDG−Pt 1STDG−Pt 2STDG−Pt 3STDG−Pt 4LSERK4tO(h )t O(h )t8O(h )tO(h )4610Fig. 2: Dissipative (up)/dispersive (down) convergence rates of theeigenvalues of the evolution operator H for different orders & ht.V. NUMERICAL RESULTSFor validation, we have simulated the above problem usinga discretization of Kx = 8, hx = 0.125, ht = 0.1, andP t = P x = 4 and upwind fluxes up to a time T = 10000(note that all quantities are dimensionless in our system ofunits). Also, we have computed the same problem with theexplicit implementation of the scheme. The STDG schemehas a single evolution operator with 6400 non-zero entries.The E-STDG works with a split domain having 8 differentevolution operators totaling 4800 non-zero entries, indicatinga significant reduction in the computational complexity. Fig.3 shows the resonances and the error for the first modes.Fig. 3: Response to a white noise in the E-STDG and STDG schemeswith P t = 4, ht = 0.1, P x = 4 ,Kx = 8, hx = 0.125 after a timeof 100. Vertical dashed lines represent the analytical modes.VI. CONCLUSIONSIn this work, we have introduced and analyzed a novelimplicit and stable nodal STDG technique, and compared itwith the classic DG-LSERK4 scheme. This nodal formulationneeds to store twice the DOFs of the LSERK4 formulations,but it is suitable to be used in high-order-in-space schemes,taking full advantage of the convergence of the spatial DGsemi-discretization. To overcome its implicitness, we haveproposed a new explicit implementation, showing that thescheme can be computationally affordable, allowing certainfreedom in choosing the size of the time step and evolutionoperators.REFERENCES[1] J. S. Hesthaven and T. Warburton, Nodal Discontinuous GalerkinMethods: Algorithms, Analysis, and Applications, 1st ed. SpringerPublishing Company, Incorporated, 2007.[2] M. Liu, K. Sirenko, and H. Bagci, “An efficient discontinuous Galerkinfinite element method for highly accurate solution of Maxwell equations,”Antennas and Propagation, IEEE Transactions on, vol. 60, no. 8, pp.3992–3998, 2012.[3] S. Petersen, C. Farhat, and R. Tezaur, “A spacetime discontinuousGalerkin method for the solution of the wave equation in thetime domain,” International Journal for Numerical Methods inEngineering, vol. 78, no. 3, pp. 275–295, 2009. [Online]. Available:http://dx.doi.org/10.1002/nme.2485[4] J. van der Vegt and H. van der Ven, “Space-time discontinuousGalerkin finite element method with dynamic grid motion for inviscidcompressible flows: I. general formulation,” Journal of ComputationalPhysics, vol. 182, no. 2, pp. 546 – 585, 2002. [Online]. Available:http://www.sciencedirect.com/science/article/pii/S0021999102971858[5] P. Monk and G. R. Richter, “A discontinuous Galerkin method forlinear symmetric hyperbolic systems in inhomogeneous media,” Journalof Scientific Computing, vol. 22-23, pp. 443–477, 2005. [Online].Available: http://dx.doi.org/10.1007/s10915-004-4132-5[6] R. S. Falk and G. R. Richter, “Explicit finite element methods forsymmetric hyperbolic equations,” SIAM J. Numer. Anal., vol. 36, no. 3,pp. 935–952, Mar. 1999. [Online]. Available: http://dx.doi.org/10.1137/S0036142997329463[7] J. C. Meza and S. W. W., “Domain decomposition algorithms for linearhyperbolic equations,” Rice University, Tech. Rep., 1987.[8] J. Alvarez, L. D. Angulo, M. R. Cabello, A. R. Bretones, andS. S. G. Garcia, “An analysis of the leap-frog discontinuous galerkinmethod for maxwell’s equations,” Microwave Theory and Techniques,IEEE Transactions on, vol. 62, no. 2, pp. 197–207, Feb 2014. [Online].Available: http://dx.doi.org/10.1109/TMTT.2013.2295775

https://digibug.ugr.es/bitstream/10481/50157/1/a4_anexplicit_mwcl14.pdf

An Explicit Nodal Space-Time Discontinuous Galerkin Method for Maxwell’s Equations

Abstract

Similar works

Full text

Available Versions

Institutional Repository of the University of Granada

Repositorio Institucional Universidad de Granada