A recently introduced time-domain method for the numerical solution of Maxwell's equations on unstructured grids is reformulated in a novel way, with the aim of implementation on graphical processing units (GPUs). Numerical tests show that the GPU implementation of the resulting scheme yields correct results, while also offering an order of magnitude in speedup and still preserving all of the main properties of the original finite-difference time-domain algorithm. © 2017 IEEE

Cicuttin, Matteo

Codecasa, Lorenzo

Kapidani, Bernard

Specogna, Ruben

Trevisan, Francesco

Archivio istituzionale della ricerca - Università degli Studi di Udine

10 April 2018Università degli studi di UdineGPU Accelerated Time-Domain Discrete Geometric Approach Method for Maxwell's Equations on Tetrahedral Grids /Cicuttin, Matteo; Codecasa, Lorenzo; Kapidani, Bernard; Specogna, Ruben; Trevisan, Francesco. - In: IEEETRANSACTIONS ON MAGNETICS. - ISSN 0018-9464. - STAMPA. - 54:3(2018), pp. 1-4.OriginalGPU Accelerated Time-Domain Discrete Geometric Approach Method for Maxwell'sEquations on Tetrahedral GridsPublisher:PublishedDOI:10.1109/TMAG.2017.2753322Terms of use:Publisher copyright(Article begins on next page)The institutional repository of the University of Udine (http://air.uniud.it) is provided by ARIC services. Theaim is to enable open access to all the world.Availability:This is the peer reviewd version of the followng article:This version is available http://hdl.handle.net/11390/1127249 since 2018-03-05T10:37:59ZIEEE TRANSACTIONS ON MAGNETICS, VOL. 00, NO. 0, JUNE 2017 1GPU accelerated time domain discrete geometric approach methodfor Maxwell’s equations on tetrahedral gridsMatteo Cicuttin1,2, Lorenzo Codecasa3, Bernard Kapidani4, Ruben Specogna4, Francesco Trevisan41Universite´ Paris-Est, CERMICS (ENPC), F-77455 Marne-la-Valle´e, France2INRIA Paris, 2 rue Simone Iff, F-75589 Paris, France3Politecnico di Milano, Dip. di Elettronica, Informazione e Bioingegneria, I-20133, Milano, Italy4Polytechnic Department of Engineering and Architecture (DPIA), Universita` di Udine, I-33100 Udine, ItalyA recently introduced time domain method for the numerical solution of Maxwell’s equations on unstructured grids is reformulatedin a novel way, with the aim of implementation on Graphical Processing Units (GPUs). Numerical tests show that the GPUimplementation of the resulting scheme yields correct results, while also offering an order of magnitude in speedup and stillpreserving all the main properties of the original Finite Difference Time Domain (FDTD) algorithm.Index Terms—Discrete Geometric Approach, Wave Propagation, Parallel computing, FDTDI. INTRODUCTIONTHE extension of Yee’s original FDTD algorithm [1] tounstructured grids is a topic of intense research in com-putational electromagnetics. Unstructured grids, differentlyfrom the cartesian grids used in the original FDTD, allow toeasily handle the complex geometries arising in applicationsfrom both experimental science and industry. However, tobe successful, FDTD-like algorithms on unstructured gridsmust retain the properties that make the original algorithmso ubiquitous. One of these is its uniform and small stencil,which makes the idea of a numerical time domain solution ofhuge wave propagation problems on machines with parallelcomputing capabilities very appealing. In practice, the algo-rithm has thus seen its performance bottleneck getting moreand more linked to the memory bandwith of the underlyinghardware architecture. For this reason the focus of FDTDcode developers has recently shifted from multicore CPUs toGraphical Processing Units (GPUs) [3], [4], which providetypically an order of magnitude more bandwidth than highperformance CPUs for the price of a medium budget laptop.On unstructured grids, a parallel implementation using acluster of GPUs has been presented in [5] based on the timedomain application of the Discontinuous Galerkin (DG) FiniteElement Method. However, DG formulations exhibit spuriousmode solutions due to violation of charge (and possiblyenergy) conservation [6].In the present paper we will use the approach of [7],based on the Discrete Geometric Approach (DGA), whichyields a conditionally stable, explicit and consistent schemeon tetrahedral grids. Furthermore the algorithm, to which wewill refer as the DGATD in the following, is stencil-basedas the original FDTD. The accuracy and performance of thisapproach on single-core CPU architectures have already beenManuscript received June ??, 2017; accepted ??. Date of current ver-sion June 27, 2017. Corresponding author: B. Kapidani (e-mail: kapi-dani.bernard@spes.uniud.it).Color versions of one or more of the figures in this paper are availableonline at http://ieeexplore.ieee.org.Digital Object Identifier 00.0000/TMAG.2017.0000000studied in [8]. The aim of this paper is instead to present areformulation of this approach which also exhibits the samepromise for exploiting hardware parallelization as the originalFDTD method.The paper is organized as follows: in Section II the spaceand time discretization approach is outlined in detail, inSection III a brief stability analysis for the resulting numer-ical scheme is given, in Section IV the correctness of theimplementation is assessed on a test case with known analyt-ical solution, in Section V implementation and performancerelated details within the NVIDIA Compute Unified DeviceArchitecture (CUDA) programming model are given. Finally,in Section VI, some conclusions are drawn.II. FRACTIONED GRID FORMULATIONWe wish to solve the initial value problem for Maxwell’sequationsµ∂h∂t= −∇× e, (1)ε∂e∂t= ∇× h, (2)in a bounded region Ω ⊂ R3, in which we assume noconductive losses. In the framework of the Discrete GeometricApproach (DGA), electromagnetic quantities are discretizedinto degrees of freedom (DoFs) by their integrals over geo-metric elements of two dual grids: a primal tetrahedral meshG and a polyhedral dual mesh (dual complex) G˜ obtained bybarycentric subdivision of G. In the DGA, discrete differentialoperators are encoded by the incidence matrices betweengeometric elements in each mesh: in (1) the flux of the curl ofe through any triangle of G equals the (oriented) sum of thecirculations of e over the edges in its boundary; similarly in(2) the flux of the curl of h through any face of G˜ equalsthe (oriented) sum of the circulations of h over the dualedges in its boundary. A remarkable consequence is that forMaxwell’s equations written in integral form the curl operatordiscretization is exact.IEEE TRANSACTIONS ON MAGNETICS, VOL. 00, NO. 0, JUNE 2017 2Fig. 1. Hexahedron Ωr is the element on which the basis functions aredefined.Let us take any disjoint intersection Ωr between a tetrahe-dron τn of G and a dual volume τ˜v of G centered in any ofthe four vertices v of τn (see Fig. 1). By construction, Ωris always a hexahedron. Furthermore, τn is the union of fourdisjoint hexahedra, while τ˜v is the union of a variable numberof disjoint hexahedra (equal to the number of tetrahedra inwhich v is a vertex). We will define fractioned DoFs in eachΩr asfˆr = h · a˜1rh · a˜2rh · a˜3r , uˆr = e · a1re · a2re · a3r ,where a˜1r , a˜2r , a˜3r are the tangent vectors of the portions of thethree edges of τ˜v in the boundary of Ωr (as oriented in Fig.1) and a1r , a2r , a3r are the tangent vectors of the portions of thethree edges of τn in the boundary of Ωr (as oriented in Fig.1). These DoFs are the main difference with respect to [7] andwill lead to a different, although equivalent, structure in theresulting algorithm. Furthermore we define a local face-edgeincidence matrix in ΩrCr = 0 −1 11 0 −1−1 1 0 ,and two distinct sets of uniform valued basis functions w˜ir(r)and wir(r), with i = 1, 2, 3, which have compact support andcan be written asw˜ir(r) =a˜jr × a˜kra˜ir × a˜jr · a˜kr, (3)wir(r) =ajr × akrair × ajr · akr, (4)in which i, j, k is any permutation of 1, 2, 3. Finally we definetwo local mass matrices, whose element at the i-th row andj-th column isMµr (i, j) =∫Ωrw˜ir(r) · µ(r)w˜jr(r) dr,Mεr (i, j) =∫Ωrwir(r) · ε(r)wjr(r) dr.Since Mµr and Mεr are 3 × 3 symmetric, positive-definitematrices we are easily tempted to discretize Eqs. (1)–(2) ineach Ωr. To impose continuity of tangential field componentsFig. 2. Procedure of assembling local mass matrices: detail of two adjacenttetrahedra.as in the original FDTD, we proceed thinking in globalvariables terms: the set of unknowns of the discrete problemcomprises the set uˆ, i.e. the circulations of e over each halfof each edge of G without repetitions, and the set fˆ , i.e. thecirculations of h over each half of each edge of G˜, againwithout repetitions (we split a full dual edge at the barycenterof the triangle it crosses). To get a consistent discretizationof (1), we have to glue together the contributions of the fourΩr-like hexahedra in τn. Take for example the simple case ofFig. 2 in which G is just two tetrahedra: the yellow and thered hexahedra are in the same tetrahedron, so their local Mµrcan be assembled to contribute to the local mass matrix of afull tetrahedron. If we do likewise for the other two Ωr-likehexahedra in the same tetrahedron we getMµτnfˆn+ 12τn − fˆn−12τn∆t= −Cτnuˆn,where Mµτn is a 4 × 4 matrix and fˆτn is a column vectorof size 4. Matrix Cτn is a (very sparse) superposition ofthe appropriate Cr matrices, and again can be constructedefficiently.Referring again to Fig. 2, the yellow and the blue hexahedraare instead in two different tetrahedra, therefore they contributeto separate Mµτn . Yet these two hexahedra share a vertexv of G, so they will contribute the a local mass matrixMετ˜v , associated to the volume of G˜ centered in v. Then,to get a consistent discretization of (2), we glue together thecontributions of all the Ωr-like hexahedra in a dual τ˜vMετ˜vuˆn+1τ˜v − uˆnτ˜v∆t= C˜τ˜v fˆn+ 12 ,where Mετ˜v is a d˜× d˜ matrix, uˆτ˜v is a column vector of sized˜, where d˜ is the number of edges of G which contain v and|v| is the number of tetrahedra of G which contain v. MatrixC˜τ˜v is a (very sparse) superposition of the appropriate setof CTr matrices and can be constructed efficiently. We notethe remarkable fact that, by requiring tangential continuityof e and h and exploting the properties of (3)–(4), matricesMµτn and Mετ˜vare locally and globally exactly divergencepreserving matrices for piecewise-constant fields e = eτ˜vand h = hτn (see [7] for the proof). Furthermore they aresymmetric positive-definite, small, sparse and can be invertedIEEE TRANSACTIONS ON MAGNETICS, VOL. 00, NO. 0, JUNE 2017 3Fig. 3. Result of the FFT on the computed electric field in the cylindricalresonator. The resonances match the theoretically predicted frequency values.locally, yieldingfˆn+ 12τn = fˆn− 12τn −∆t(Mµτn)−1Cτnuˆn ∀τn ∈ G, (5)uˆn+1τ˜v = uˆnτ˜v + ∆t(Mετ˜v )−1C˜τ˜v fˆn+ 12 ∀τ˜v ∈ G˜, (6)which is a valid, explicit, time marching scheme. We remarkthat the local curl matrices Cτn , C˜τ˜v are very sparse and theirproduct with global field vectors can be performed in parallel.We will refer to (5)–(6) as the fractioned Discrete GeometricApproach Time Domain (fractioned DGATD) method in therest of the paper. We remark that (5)–(6) is a reformulationof the scheme of [7], and that the two yield exactly the samesolution up to machine precision.III. STABILITY ANALYSISThe fractioned DGATD is conditionally stable due to theproperties of its mass matrices [7]. Simulations with more than300 000 time steps do not show any late time instabilities, asalso the test case of section IV will show. If we pack thesystem of equations in (5)–(6) by appending all the local massmatrices into global block-diagonal, matrices Mµ−1τ , Mε−1τ˜ ,by appending all the local curl matrices, and by appending allthe field unknowns in column vector format, the scheme readsfˆ n+12 = fˆ n−12 −∆tMµ−1τ Cτ uˆn, (7)uˆn+1 = uˆn + ∆tMε−1τ˜ C˜τ˜ fˆn+ 12 , (8)where we have dropped the subscripts on τ and τ˜ to stressthat the constitutive matrices are not restricted any more to aparticular primal or dual volume, respectively. Likewise Cτand C˜τ˜ are obtained by appending all the Cτn and Cτ˜v ,respectively. Eq. (8) can be rewritten asuˆn+1 = 2uˆn − uˆn−1 + ∆tMε−1τ˜ C˜τ˜(fˆ n+12 − fˆ n− 12),from which, using (7), it ensuesuˆn+1 =(2I−∆t2Mε−1τ˜ C˜τ˜Mµ−1τ Cτ)uˆn − uˆn−1.The timestepping structure is thus shown to be equivalent tothe one proposed on standard Finite Elements in [9]. For agiven grid, the fractioned DGATD scheme is stable providedthat ∆t < 2/√ρK, where ρK is the spectral radius of the ma-trix K = Mε−1τ˜ C˜τ˜Mµ−1τ Cτ . This means that the maximumFig. 4. Zoom in frequency. The theoretically predicted peak at 245.06 MHzis below the noise floor if the simulation is too short.allowed ∆t can be efficiently estimated using algorithms basedon the power iteration method.IV. CORRECTNESS ANALYSISThe proposed formulation was tested on a problem withknown analytical solution inspired from [10]. We tested thealgorithm on a cylindrical resonator (height h = 0.5m andradius r = 1m) filled with air. The resonator is excited viaexternal coupling with a broad Gaussian pulse. Thus we wishto solve an eigenvalue problem in the time domain, since it isknown that the resonant frequencies are given byfmnp =c2pi√µrεr√(χmnr)2+(ppih)2,where χmn is the m-th zero of the n-th cylindrical Besselfunction, p is a non-negative integer and c is the speed oflight in vacuum. We can retrieve spectral behaviour (Figs. 3,4)for an arbitrarily large frequency interval after a single timedomain simulation by computing the FFT on the interpolated ein one element of G. We know from the analysis that we shouldsee spectral peaks in Fig. 3 at frequencies f1 = 114.75 MHz,f2 = 182.84 MHz, f3 = 245.06 MHz, f4 = 263.38 MHz. Thepeaks at f1, f2, and f4 are visible in Fig. 3. The peak at f3seems to be missing, but is actually retrieved by zooming in therange between 240 and 270 MHz (see Fig. 4). We remark that,to resolve the peak at f3 from the nearest peak, a rather highsampling frequency is needed, even after applying Hammingwindowing on the input waveform. This constraint translatesinto a longer simulation (a further motivation for exploitingparallelism).V. PERFORMANCE ANALYSISThe method was implemented in C++. For the single-coreCPU version we used the original DGATD of [7], running on aXeon E5-2687Wv4 processor with Eigen 3.3.1 for linear alge-bra operations. The CUDA C++ code was tested on a TESLAC2075 accelerator. Two reasonable approaches are possiblein implementing the algorithm on the GPU. The first, moreabstract approach, is to still use the original formulation of[7], which features huge sparse matrix vector multiplicationsas its only conspicuous floating point operation, and adapt it touse the cuSparse library provided by NVIDIA, which providesa sparse matrix vector multiplication (SpMV) function. On theIEEE TRANSACTIONS ON MAGNETICS, VOL. 00, NO. 0, JUNE 2017 4 0.001 0.01 0.1 1854752 2625360C omp ut at i on  t i me  ( s )Degrees of freedomCPUGPU SpMvGPU CustomFig. 5. Average computational cost of a single time step vs. number of DoFsof the problem. x axis in linear scale, y axis in logarithmic scale.test case described in Section IV, somewhat surprisingly, thisrather naive GPU implementation already provides an orderof magnitude of speedup with respect to CPU based one,as shown in of Table I. However, SpMV produces memoryaccess patterns that are not optimal on GPUs [11]. Thisdetail was actually the motivation for the reformulation ofthe algorithm presented in Section II, in which we try toexploit the block-diagonal structure of the mass matrices toimprove performance. The key insight is the following: highlatency access to the global memory of the GPU is optimizedif consecutive threads access consecutive locations of globalmemory. By looking at the right-hand side of Eq. (5) weare presented with a set of local 4×6 matrices which canbe stored using column-major ordering. This allows adjacentmatrix values to be loaded in shared memory by adjacentthreads, obtaining full memory coalescing. Local matrix vectorproducts are then computed using one thread per row.TABLE ISPEEDUPSDoFs CPU SpMV Cust854752 1 13.42 16.582625360 1 15.84 19.50Experiments on our accelerator, which provides a memorybandwidth of roughly 102 GB/s, show that a matrix multiplica-tion kernel based on this insight reaches a sustained rate of 98GB/s. The global speedup achieved with this implementationis shown again in Table I. Unfortunately the same approachdoes not perform equally well on Eq. (6), as the variable sizeof the local mass matrix makes the logic to be added trickierand heavier. For the aims of this paper, we decided to settlefor an implementation exploiting a specific kernel only forEq. (5), relying on cuSparse for the remaining computations.In Fig. 5 we show the average single time step computationaltime versus the number of unknowns of the problem for twodifferent meshes. In Fig. 6 we show the speedup achieved byusing the fractioned grid formulation in Eq. (5), with respectto the original formulation of [7]. We remark that the fact thatthere is a speedup is non-trivial, since by using the fractionedDGATD formulation we are somewhat increasing the numberof DoFs in the problem. 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014854752 2625360C omp ut at i on  t i me  ( s )Degrees of freedomGPU SpMvGPU CustomFig. 6. Average computational cost of a single execution of the kernel forFaraday’s law vs. number of DoFs of the problem. Both axes in linear scale.VI. CONCLUSIONSWe presented a reformulation of the DGATD algorithm ontetrahedral grids. The formulation is amenable to paralleliza-tion on a GPU, achieving an order of magnitude of speedup.We are currently investigating improved data layouts to furtherincrease performance; those improvements will be the subjectof a subsequent work.ACKNOWLEDGEMENTThis work was partially supported by the Strategic Plan ofthe Polytechnic Department of Engineering and Architecture(DPIA) of University of Udine.REFERENCES[1] K.S. Yee, Numerical solution of initial boundary value problems involvingMaxwells equations in isotropic media, IEEE Transactions on Antennasand Propagation, vol. AP-14, no. 3, pp. 302307, May 1966.[2] A. Taflove and S. Hagness, Computational Electromagnetics, The FiniteDifference Time Domain Method, Second Edition. Boston, MA: ArtechHouse, 2000.[3] P. Sypek, A. Dziekonski, and M. Mrozowski, How to Render FDTDComputations More Effective Using a Graphics Accelerator, IEEETransactions on Magnetics, vol. 45, no. 3, pp. 1324-1327, 2009.[4] D. De Donno, A. Esposito, L. Tarricone, L. Catarinucci, Introductionto GPU Computing and CUDA Programming: A Case Study on FDTD,IEEE Antennas and Propagation Magazine, Vol. 52, No. 3, pp. 116-122,2010.[5] N. Goedel, N. Gunn, T. Warburton, M. Clemens, Scalability ofHigher-Order Discontinuous Galerkin FEM Computations for SolvingElectromagnetic Wave Propagation Problems on GPU Clusters, IEEETransactions on Magnetics, Vol. 46, No. 8, pp. 3469-3472, 2010.[6] E. Gjonai, T. Lau, T. Weiland, Conservation Properties of the Discontin-uous Galerkin Method for Maxwell Equations, International Conferenceon Electromagnetics in Advanced Applications (ICEAA), Torino, Italy,17-21 September 2007.[7] L. Codecasa, M. Politi, Explicit, Consistent and Conditionally StableExtension of FD-TD by FIT, IEEE Transactions on Magnetics, Vol. 44,pp. 1258-1261, 2008.[8] B. Kapidani, L. Codecasa, M. Cicuttin, R. Specogna, F. Trevisan, Acomparative performance analysis of time-domain formulations for wavepropagation problems, 2016 IEEE Conference on Electromagnetic FieldComputation (CEFC), Miami, Fl., 13-16 November 2016.[9] J.M. Jin, The Finite Element Method in Electromagnetics, Second Edition,Wiley, New York, 2002.[10] M. Cinalli, A. Schiavoni, A Stable and Consistent Generalizationof the FDTD Technique to Nonorthogonal Unstructured Grids, IEEETransactions on Antennas and Propagation, Vol. 54, No. 5, pp. 1503-1512, 2006.[11] N. Bell, M. Garland, Efficient Sparse Matrix-Vector Multiplication onCUDA, NVIDIA Technical Report NVR-2008-004, Dec. 2008

GPU Accelerated Time-Domain Discrete Geometric Approach Method for Maxwell's Equations on Tetrahedral Grids

Institutional Research Information System

GPU Accelerated Time-Domain Discrete Geometric Approach Method for Maxwell&apos;s Equations on Tetrahedral Grids

A recently introduced time-domain method for the numerical solution of Maxwell's equations on unstructured grids is reformulated in a novel way, with the aim of implementation on graphical processing units (GPUs). Numerical tests show that the GPU implementation of the resulting scheme yields correct results, while also offering an order of magnitude in speedup and still preserving all of the main properties of the original finite-difference time-domain algorithm

Archivio istituzionale della ricerca - Politecnico di Milano

https://air.uniud.it/bitstream/11390/1127249/1/manuscript_dgagpu_td2017.pdf

GPU Accelerated Time-Domain Discrete Geometric Approach Method for Maxwell's Equations on Tetrahedral Grids

Abstract

Similar works

Full text

Available Versions

Archivio istituzionale della ricerca - Università degli Studi di Udine

Institutional Research Information System

Archivio istituzionale della ricerca - Politecnico di Milano