Perfectly matched layer (PML) [1, 2] is often used in finite difference , finite element , and discontinuous Galerkin time-domain (DGTD) methods [5, 6, 7, 8, 9, 10] to imitate/approximate the radiation boundary condition (i.e., truncate a unbounded physical domain to a finite computation domain) while solving Maxwell equations or the wave equation.
The performance of the PML depends on the attenuation coefficient (which is implemented as conductivity in Maxwell equations) and the thickness of the layer. Increasing either one or both of them increases the absorption inside the PML. However, in practice, one cannot use a constant and high conductivity as it would increase the numerical reflection at the interface between the PML and computation domain, or use a very thick layer since it would increase the computational cost. Therefore, a smoothly increasing conductivity profile is often used to achieve both high absorption and small numerical reflection [1, 2, 11, 3, 4, 12, 13].
In DGTD, the PML conductivity profile can be implemented in two different ways. The first method assumes that the conductivity is constant in a given element. This constant might for example be set to the value of the conductivity profile at the center of the element. Implementation of this method is rather straightforward since the mass matrices of different elements only differ by a constant (for linear elements) . However, the conductivity profile becomes discontinuous between neighboring elements [Fig. 1 (a)] and the element surfaces that are not parallel to the PML interface lead to large reflections and destroy the high-order accuracy of the solution. One workaround is to build layered (tetrahedral) meshes [Fig. 1 (b)] or use orthogonal (hybrid) meshes inside the PML [7, 14], and accordingly set a layered conductivity profile. But this makes the setups of the computation domain and the PML rather tedious since one needs to control the mesh and the conductivity on all face, edge, and corner regions of the PML. Moreover, to reduce the numerical reflection by decreasing the conductivity discontinuity between neighboring mesh (or conductivity) layers, their number has to be increased.
The second method allows the conductivity to vary inside a given element (following the increasing conductivity profile inside the PML). It should be noted here that (higher-order) DGTD allows for sampling of the material properties at the sub-elemental level. In this case, the behavior of the PML is determined only by the conductivity samples within the elements and therefore the mesh interfaces can be aligned arbitrarily [Fig. 1 (c)] without adversely affecting its performance. Indeed, this second approach can improve the PML performance [15, 16, 17]. However, the conductivity varying within the elements results in an element-dependent mass matrix. Therefore, a direct implementation requires a different mass matrix (or its inverse) to be stored for every element [9, 7, 18, 15, 16, 14]. This significantly increases DGTD’s memory footprint. For example, for the stretched-coordinate (SC)-PML , the memory required to store the mass matrices scales with , where
is the number of interpolating nodes in each element,is the number of elements in the PML, and
comes from the five material-dependent coefficients in the update equations and three Cartesian components of the vector field. In contrast, for the first method, where the conductivity is assumed constant in a given element, this memory requirement scales withsince only the constant conductivity of each element and a single reference mass matrix are stored.
In this work, a memory-efficient method to implement the SC-PML with smoothly-varying attenuation coefficients in DGTD is developed. The proposed method allows the conductivity to vary inside the elements and constructs the resulting local mass matrices using a weight-adjusted approximation (WAA) . Compared to the direct implementation that is briefly described above, the proposed method reduces the memory requirement to , where , while maintaining the PML performance. It should be noted here that the WAA has been proved to be energy-stable and preserve the high-order convergence of DGTD [19, 20, 21]. Indeed, numerical examples presented here also show the proposed method maintains the higher-order accuracy of the solution. Additionally, it is demonstrated that the PML with smoothly-increasing conductivity profile as implemented with the proposed method performs better than the PML implemented using element-wise constant conductivity profile.
Ii-a WAA-DGTD for SC-PML
Maxwell equations in stretched-coordinates for a source-free and lossless medium can be expressed as 
where and are the electric and magnetic fields, and are the permittivity and the permeability, is the frequency, and
where and are one-dimensional positive real functions along direction . Here, is the attenuation coefficient that ensures absorption inside the SC-PML and changes the propagation speed (and also increases the absorption for evanescent waves). It is well known that using smoothly-increasing and reduces the numerical reflection from the interface between the PML and the computation domain while maintaining a high absorption rate, and therefore improves the PML performance [11, 12, 13].
The time-domain update equations in SC-CPML can be expressed as 
are diagonal tensors with entries defined as
In (9) and the rest of the text follows the permutation .
Following the standard nodal discontinuous Galerkin method , the computation domain and the PML are discretized into elements with volumetric support and boundary surface , and in each element , , and are expanded using the Lagrange polynomials , , where is the number of interpolating nodes and is the order of the Lagrange polynomials. Finally, Galerkin testing yields the semi-discrete system of equations as
Here, , , , and are vectors storing the unknown time-dependent coefficients of the relevant basis functions, and , , are the mass matrices with entries
denotes the curl operator with its component along direction is given by
where , , is the component of the numerical flux along direction , which in general involves unknowns from the current element and its neighboring element [8, 22, 23], and are the stiffness and the face mass matrices with entries
and and are the permittivity and the permeability (assumed constant) in each element. Note that, the nodal DG framework  is used here, but the proposed method can be easily extended to vector DG methods [6, 7, 9, 24, 25, 26].
For linear elements, the mass matrices in (14) are simply scaled versions of the mass matrix defined on the reference element, , where is the Jacobian of the coordinate transformation between element and the reference element. Hence, only and (scalar constant) are stored. Similarly, in (15), if is assumed constant inside the elements, i.e., , then and one only needs to store in addition to and . In this case, (10)-(13) can be implemented as efficiently as the case without the PML [8, 27, 28].
However, if is allowed to vary inside the elements, are different in different elements and in general there is no simple relationship between these different mass matrices. One has to store every one of these mass matrices (or their inverse). As an alternative, the mass matrix can be recomputed at each time step, but this would significantly increase the cost of time marching . The memory required to store the mass matrices in (10)-(13) scales with per element, where is the number of the components of the vector field, is the number of the coefficients , , , , and . Note that this memory requirement is significantly higher than that of storing the unknowns coefficients of the basis functions, which scales with in the PML.
To reduce the memory requirement of implementing (15) with allowed to vary inside the elements, WAA  is used. It has been shown that with this approximation DGTD retains provable energy-stability and high-order accuracy [19, 20, 21]. Note that in the above SC-PML formulation, directly multiplying (5) and (6) with on both sides reduces the number of element-dependent mass matrices to . But this would result in a non-conservative form, whose solution is neither provably energy-stable nor provably high-order accurate [8, 19].
First, a weight-adjusted inner product is introduced to approximate the parameter-weighted inner product in the expression of the mass matrix . The mass matrix, which is associated with the element and a locally varying coefficient , is approximated as
where , , are the Gaussian quadrature nodes corresponding to the quadrature rules of degree and are the corresponding weights. Hence,
where is an interpolation matrix defined on the reference element, , and are generalized Vandermonde matrices with entries and , respectively, is the -th orthonormal polynomial basis , is also element-independent, and is a diagonal matrix containing the coefficients evaluated at the quadrature nodes.
Here, is introduced to simplify the implementation. In (21), and are defined on the reference element, and is a scaled version of the reference matrix, .
Ii-B Computational complexity
In DGTD with explicit time marching, all operations are localized within the elements. The memory required to store the mass matrices in the direct implementation of (10)-(13) scales with , where comes from the number of unknown components times the number of different mass matrices associated with different coefficients and is the number of elements in the PML. In the WAA formulation (26)-(29), the memory requirement reduces to , where comes from the number of unknown components times the number of coefficient samples at the quadrature points and comes from and defined on the reference element. For simplicial quadrature rules that are exact for up to polynomials of degree , [31, 32].
To compare the number of arithmetic operations required by the two implementations, one should first note that the curl operator is the same in both formulations. Computation of requires those of the spatial derivatives and the numerical flux [33, 34, 35]. Here, the memory access time is much more significant than the time required to carry out these computations because data from neighboring elements, which are discontinuous in memory, is required. Therefore, only the times required to complete the arithmetic operations of the remaining terms are compared. For the same reason, in practice, the time required to compute
dominates the overall time required by the time marching, and the difference in the numbers of arithmetic operations as estimated below for the remaining terms is less significant (see the example in SectionIII).
In (10), the three matrix-vector multiplications and two vector-vector additions require multiplication operations and addition operations, respectively. In (26), the multiplication of with a vector of length , and the multiplication of with a vector of length require multiplication operations. The multiplication of a diagonal matrix with a vector (such as ) requires multiplication operations. As a result, excluding the computation of , (26) requires multiplications and additions. For the auxiliary variable, the cost of (12) is multiplications and subtractions, while (28) requires multiplications and subtractions. One can see the number of operations in the WAA implementation is slightly higher than that in the direct implementation. But as mentioned above the time required by these operations is smaller than the time required to compute , and therefore overall times required by the two implementations are not that different.
Iii Numerical Examples
In this section, the accuracy and the efficiency of the proposed WAA formulation are compared to those of the traditional PML implementations using numerical examples. To this end, four PML configurations/implementations are considered in these examples: (i) and/or , , are assumed constant inside the elements on a paved mesh (EC-paved) [Fig. 1(a)], (ii) and/or , , are assumed constant inside the elements on a layered mesh (EC-layered) [Fig. 1(b)], (iii) and/or , , are allowed to vary inside the elements on a paved mesh (SV-paved) [Fig. 1(c)], and (iv) same configuration in (iii) but implemented using the proposed method with the WAA (SV-WAA-paved) [Fig. 1(c)]. In all implementations, the order of the Lagrange polynomials , which results in .
For configuration (i), the constant values in a given element are obtained by sampling and , , at that element’s node that is farthest away from the PML interface (along the -direction). For configuration (ii), to ensure that the element surfaces are strictly parallel to the axes, the PML mesh is built layer by layer and constant values in a given layer are obtained by sampling and , , at the outermost surface of that layer (along the -direction). For the WAA in implementation (iv), the order of the Gaussian quadrature rule is , resulting in .
In all examples, the background medium is free space and the excitation is a plane wave with electric field , where , is the speed of light in free space, and is a base-band Gaussian pulse with and . The average edge lengths of all meshes used under this excitation are .
First, the reflection of a plane wave normally incident on the PML is computed. The computation domain is a rectangular box with dimensions . Perfect electric conductor (PEC) and periodic boundary conditions are used on the outer boundary of the PML that is located perpendicular to the direction and on the computation domain boundaries perpendicular to the and directions, respectively. The plane wave excitation is introduced on surface and propagates in the -direction. The domain is long enough to ensure that the reflected field is well-separated from the incident one, and therefore the reflection from the PML is simply measured by the peak value of the reflected field’s amplitude. The conductivity profile is described by , where is the -coordinate on the interface between PML and the computation domain, is the thickness of the PML and is the order of the profile. Note that is nonzero only when . The values of these parameters are , , and , and also both inside the PML and the computation domain.
In this example, four configurations/implementations are considered: EC-paved, EC-layered, SV-paved, and SV-WAA-paved. Their performances are compared for . For all four groups of simulations, is scanned to find the minimum reflection that can be obtained for each case. Fig. 2 shows that with increasing , the reflection first decreases exponentially and then increases gradually. This is observed for all configurations/implementations and all values of . When is small, the overall reflection is dominated by the reflection from the PEC boundary simply because the absorption inside the PML is not high enough. Therefore, in this regime, increasing elevates the absorption and reduces the amplitude of the wave reflected back into the computation domain exponentially. The numerical reflection (which is smaller than the reflection from the PEC boundary for small ) increases with increasing , and starts dominating the overall reflection as demonstrated in the figure by the gradual increase after the minimum point.
For the EC-paved configuration, the reflection stays at a high level and does not decrease with increasing . This is because of the large reflections from unoriented internal element surfaces. For the EC-layered, SV-layered, and SV-WAA-paved configurations, high-order convergence is observed, i.e., the reflection keeps on decreasing exponentially with increasing . Still, the reflection for the SV-paved and SV-WAA-paved configurations is about smaller than that for the EC-layered configuration. Note that this higher accuracy comes with the ease of meshing since a layered mesh (and conductivity profile) is not needed. Finally, Fig. 2 also shows that the SV-WAA-paved implementation performs exactly the same as the SV-paved direct implementation, which verifies the accuracy of the proposed method.
Next, scattering from a PEC sphere of radius is considered. The computation domains and the PMLs for the EC-layered and SV-paved and SV-WAA-paved configurations are shown in Figs. 3 (a) and (b), respectively. The plane wave excitation is introduced on the total-field scattered-field (TFSF) surface [shown in green in Figs. 3 (a) and (b)]. The conductivity function is , , , where is the -coordinate on the interface between PML and the computation domain, is the thickness of the PML along the direction, and is the order of the profile. Note that is nonzero only when . The values of these parameters are , and .
Because the distance between the sphere surface and the PML is short, possibly-evanescent scattered waves enter the PML with high grazing angles. A varying , , profile is employed to help with the absorption of these evanescent waves [12, 13]: with . Note that inside the computation domain, . In this example, using the PEC or the first-order absorbing boundary condition  on the outer boundary of the PML gives similar results. The results presented here are obtained with the PEC boundary condition.
Three configurations/implementations are considered here: EC-layered, SV-paved, and SV-WAA-paved. For the EC-layered configuration, the thickness of each PML mesh layer is [Fig. 3 (a)]. Note that, for this example, generation of these layers is rather tedious since in the corner region one has to align all layer/element surfaces in all three directions. In contrast, for the SV-paved and SV-WAA-paved configurations [same mesh is used – Fig. 3 (b)], and values are simply obtained by sampling the corresponding profile functions at the nodes of the elements. This significantly simplifies the setups of the computation domain and the PML since even an explicit interface between the computation domain and the PML is not required [see Fig. 3 (b)]. The performances of the three configurations are compared for and .
is scanned to find the minimum reflection that could be reached for each case. Note that in this example “reflection” is defined as the peak value of the absolute difference between the fields computed at a probe point for the above cases and corresponding reference fields computed at the same point. different probe points (placed either in the TF or in the SF region) have been tested and the results are consistent for all of them. The results below correspond to the probe point at . These reference fields are computed under the same excitation but the distance between the sphere surface and the PML is extended to . To ensure that the discretization errors are at the same level, the meshes in the overlapped regions between the actual computation domains and the extended ones are kept exactly the same, the average edge lengths of the meshes in the extended region are kept the same as those in the actual computation domains, and the solutions are obtained using the same value of .
Fig. 4 plots the reflection for the three cases with different values of and versus . Clearly, the SV-paved and SV-WAA-paved configurations perform better than the EC-layered configuration for every value of . The best performance is obtained with . Note that further increasing degrades the PML performance for all configurations/implementations and all values of since high conductivity values only appear at the very end of the PML when is high. Fig. 4 also shows that the SV-WAA-paved implementation performs exactly the same as the SV-paved direct implementation, which means the error caused by the WAA of the mass matrices is below the level of the discretization error.
|memory (KB)||CPU time per step (s)|
Tested on a workstation with Intel Xeon(R) E5-2680 v4 CPU and 128GB memory. A single process is used. and .
Table. I compares the computational cost of the SV-paved and SV-WAA-paved implementations. With increasing , the memory requirement increases dramatically for the SV-paved direct implementation but only modestly for the SV-WAA implementation. For , the memory requirement of the SV-paved implementation is times that of the SV-WAA-paved implementation. The computation time required by the SV-WAA-paved implementation per time step is slightly larger than that required by the SV-paved implementation due to the increased number of arithmetic operations (see Section. II-B). It should be also be noted here that, in practice, a DGTD algorithm is usually parallelized. The difference in times required for updating different elements is relatively small and can be easily compensated by allocating a smaller number of elements for those MPI processes containing PML elements. In the numerical results presented here, assigning a weight of for PML elements in ParMetis [36, 37] yields a good load-balance.
A PML implementation that allows the attenuation coefficient to vary inside the discretization elements yields a smaller numerical reflection from the interface between the PML and the computation domain and significantly simplifies the meshing process. However, these advantages come at the cost of increased memory footprint since a different mass matrix has to be stored for every discretization element. In this work, this memory requirement is reduced by applying WAA to the mass matrices without abandoning the advantages listed above. Indeed, numerical results demonstrate that the PML with smoothly-increasing conductivity profile as implemented with the proposed method performs better than the PML implemented using element-wise constant conductivity profile and that the higher-order accuracy of the solution is maintained.
The proposed method is especially useful for simulations running on shared-memory systems where the high memory requirement of smoothly-varying PMLs could be a bottleneck. For simulations running on distributed-memory systems, the memory requirement of a single computing node is also reduced and a better load-balance could be reached with a slightly adjusted weight in the domain partition.
-  J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys., vol. 114, no. 2, pp. 185–200, 1994.
-  W. C. Chew and W. H. Weedon, “A 3D perfectly matched medium from modified Maxwell’s equations with stretched coordinates,” Microw. Opt. Tech. Lett., vol. 7, no. 13, pp. 599–604, 1994.
-  A. Taflove and S. C. Hagness, Computational electrodynamics: the finite-difference time-domain method. Artech house, 2005.
-  J.-M. Jin, The finite element method in electromagnetics. John Wiley & Sons, 2015.
-  J. Hesthaven and T. Warburton, “Nodal high-order methods on unstructured grids: I. time-domain solution of Maxwell’s equations,” J. Comput. Phys., vol. 181, no. 1, pp. 186 – 221, 2002.
-  B. Cockburn, F. Li, and C.-W. Shu, “Locally divergence-free discontinuous Galerkin methods for the Maxwell equations,” J. Comput. Phys., vol. 194, no. 2, pp. 588 – 610, 2004.
-  T. Lu, P. Zhang, and W. Cai, “Discontinuous Galerkin methods for dispersive and lossy Maxwell’s equations and PML boundary conditions,” J. Comput. Phys., vol. 200, no. 2, pp. 549–580, 2004.
-  J. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. NY, USA: Springer, 2008.
-  S. D. Gedney, C. Luo, J. A. Roden, R. D. Crawford, B. Guernsey, J. A. Miller, T. Kramer, and E. W. Lucas, “The discontinuous Galerkin finite-element time-domain method solution of Maxwell’s equations,” Appl. Comput. Electromagn. Soc. J., vol. 24, no. 2, p. 129, 2009.
-  G. C. Cohen and S. Pernet, Finite element and discontinuous Galerkin methods for transient wave equations. Springer, 2017.
-  W. C. Chew and J. Jin, “Perfectly matched layers in the discretized space: An analysis and optimization,” Electromagn., vol. 16, no. 4, pp. 325–340, 1996.
-  J.-P. Berenger, “Perfectly matched layer (PML) for computational electromagnetics,” Synthesis Lectures on Computational Electromagnetics, vol. 2, no. 1, pp. 1–117, 2007.
-  S. D. Gedney, “Introduction to the finite-difference time-domain (fdtd) method for electromagnetics,” Synthesis Lectures on Computational Electromagnetics, vol. 6, no. 1, pp. 1–250, 2011.
-  G. Chen, L. Zhao, W. Yu, S. Yan, K. Zhang, and J.-M. Jin, “A general scheme for the discontinuous Galerkin time-domain modeling and s-parameter extraction of inhomogeneous waveports,” IEEE Trans. Microwave Theory Tech., vol. 66, no. 4, pp. 1701–1712, 2018.
-  L. Angulo, J. Alvarez, M. Pantoja, S. Garcia, and A. Bretones, “Discontinuous Galerkin time domain methods in computational electrodynamics: State of the art,” in Forum Electromagn. Res. Methods Appl. Technol., vol. 10, 2015, pp. 1–24.
-  L. M. D. Angulo, “Time domain discontinuous Galerkin methods for Maxwell equations,” Ph.D. dissertation, Universidad de Granada, 2014.
-  K. Sankaran, “Accurate domain truncation techniques for time-domain conformal methods,” Ph.D. dissertation, ETH Zurich, 2007.
-  J. Niegemann, M. Konig, K. Stannigel, and K. Busch, “Higher-order time-domain methods for the analysis of nano-photonic systems,” Photonic. Nanostruct., vol. 7, no. 1, pp. 2–11, 2009.
-  J. Chan, R. J. Hewett, and T. Warburton, “Weight-adjusted discontinuous Galerkin methods: wave propagation in heterogeneous media,” SIAM J. Sci. Comput., vol. 39, no. 6, pp. A2935–A2961, 2017.
-  K. Guo and J. Chan, “Bernstein-Bezier weight-adjusted discontinuous Galerkin methods for wave propagation in heterogeneous media,” J. Comput. Phys., vol. 400, p. 108971, 2020.
-  K. Shukla, J. Chan, V. Maarten, and P. Jaiswal, “A weight-adjusted discontinuous galerkin method for the poroelastic wave equation: penalty fluxes and micro-heterogeneities,” J. Comput. Phys., vol. 403, p. 109061, 2020.
-  L. Chen and H. Bagci, “Steady-state simulation of semiconductor devices using discontinuous Galerkin methods,” IEEE Access, vol. 8, pp. 16 203–16 215, 2020.
-  L. Chen and H. Bagci, “Multiphysics modeling of plasmonic photoconductive devices using discontinuous Galerkin methods,” arXiv preprint arXiv:1912.03639, 2019.
-  P. Li, Y. Shi, L. J. Jiang, and H. Bagci, “DGTD analysis of electromagnetic scattering from penetrable conductive objects with IBC,” IEEE Trans. Antennas Propag., vol. 63, no. 12, pp. 5686–5697, 2015.
-  P. Li, L. J. Jiang, and H. Bagci, “Transient analysis of dispersive power-ground plate pairs with arbitrarily shaped antipads by the DGTD method with wave port excitation,” IEEE Trans. Electromagn. Compat., vol. 59, no. 1, pp. 172–183, 2017.
-  ——, “Discontinuous Galerkin time-domain modeling of graphene nanoribbon incorporating the spatial dispersion effects,” IEEE Trans. Antennas Propag., vol. 66, no. 7, pp. 3590–3598, 2018.
-  M. Liu, K. Sirenko, and H. Bagci, “An efficient discontinuous Galerkin finite element method for highly accurate solution of Maxwell equations,” IEEE Trans. Antennas Propag., vol. 60, no. 8, pp. 3992–3998, 2012.
-  L. Chen and H. Bagci, “A discontinuous Galerkin framework for multiphysics simulation of photoconductive devices,” in Proc. Int. Appl. Comput. Electromagn. Symp. IEEE, 2019, pp. 1–2.
-  K. Sirenko, M. Liu, and H. Bagci, “Incorporation of exact boundary conditions into a discontinuous Galerkin finite element method for accurately solving 2d time-dependent Maxwell equations,” IEEE Trans. Antennas Propag., vol. 61, no. 1, pp. 472–477, 2012.
-  L. Chen and H. Bagci, “A unit-cell discontinuous Galerkin scheme for analyzing plasmonic photomixers,” in Proc. IEEE Int. Symp. Antennas Propag., 2019, pp. 1069–1070.
-  R. Cools, “Monomial cubature rules since “stroud”: A compilation – part 2,” J. Comput. Appl. Math., vol. 112, no. 1-2, pp. 21–27, Nov. 1999.
-  H. Xiao and Z. Gimbutas, “A numerical algorithm for the construction of efficient quadrature rules in two and higher dimensions,” Comput. Math. with Appl., vol. 59, no. 2, pp. 663–676, 2010.
-  L. Chen, M. Dong, and H. Bagci, “Modeling floating potential conductors using discontinuous Galerkin method,” IEEE Access, vol. 8, pp. 7531–7538, 2020.
-  K. Sirenko, Y. Sirenko, and H. Bagci, “Exact absorbing boundary conditions for periodic three-dimensional structures: Derivation and implementation in discontinuous Galerkin time-domain method,” IEEE J. Multiscale and Multiphys. Comput. Techn., vol. 3, pp. 108–120, 2018.
-  L. Chen, M. Dong, P. Li, and H. Bagci, “A hybridizable discontinuous Galerkin method for simulation of electrostatic problems with floating potential conductors,” arXiv preprint arXiv:2006.02549, 2020.
-  G. Karypis and V. Kumar, “A parallel algorithm for multilevel graph partitioning and sparse matrix ordering,” J. Parallel Distr. Com., vol. 48, no. 1, pp. 71–95, 1998.
-  L. Chen and H. Bagci, “An MPI-based parallel multiphysics discontinuous Galerkin framework for photoconductive devices,” in Prog. in Electromagn. Res. Symp. December 17-20, 2019 in Xiamen, China, 2019. [Online]. Available: http://hdl.handle.net/10754/662548