Efficient Time Integration of Maxwell's Equations with Generalized Finite Differences

. We consider the computationally eﬃcient time integration of Maxwell’s equations using discrete exterior calculus (DEC) as the computational framework. With the theory of DEC, we associate the degrees of freedom of the electric and magnetic ﬁelds with primal and dual mesh structures, respectively. We concentrate on mesh constructions that imitate the geometry of the close packing in crystal lattices that is typical of elemental metals and intermetallic compounds. This class of computational grids has not been used previously in electromagnetics. For the simulation of wave propagation driven by time-harmonic source terms, we provide an optimized Hodge operator and a novel time discretization scheme with nonuniform time step size. The numerical experiments show a signiﬁcant improvement in accuracy and a decrease in computing time compared to simulations with well-known variants of the ﬁnite diﬀerence time domain method.

1. Introduction.We consider the time-dependent Maxwell's system formulated in a three-dimensional exterior space domain Ω as (E, B, H, D) (t 0 ) = E 0 , B 0 , H 0 , D 0 , in Ω, (1.3) where E and H are the electric and magnetic fields, and D and B are the electric and magnetic flux densities.The initial conditions, E 0 , B 0 , H 0 , and D 0 are set at the initial time t 0 , while the final state is considered at the final time t tot .The system is applied with the constitutive relations D = E and B = µH, where electric permittivity, , and magnetic permeability, µ are tensor-valued space-dependent variables.The source functions, J = −σE and J * = −σ * H, depend on the tensor-valued space-dependent variables electric conductivity, σ, and magnetic conductivity, σ * .
For solving (1.1)-(1.3)numerically, the physical quantities and the differential operators must be discretized.Conventional methods, such as the standard finite difference or finite element method, pay no attention to the geometric structure of the physics behind the models.Since the vector-valued electric and magnetic fields are not well-presented by nodal values, non-physical modes may occur in the numerical solution [28].That is why there is a growing interest in methods that maintain the fundamental properties of continuum equations at the discrete level [1].
The finite difference time domain (FDTD) method, introduced by Yee [62], has been used for several decades as a common method for solving electromagnetic problems in the space-time domain.Later on, it was recognized that the method's fashion of locating the electric field variables on the edges and the magnetic field variables at the centers of the mesh element surfaces is the key concept for preserving the continuity of physical phenomena [25].Originally, the FDTD was restricted to the structured grids constructed by cubic elements, but there are certain generalizations to, e.g., three-dimensional nonuniform [39] and non-orthogonal [32] grids and two-dimensional polyhedra [27].Other physics-preserving methods that work in the style of staggered finite difference grids include, e.g., mimetic finite difference (MFD) schemes [34].
The above-mentioned mentioned generalizations of the FDTD can be covered by employing the metric-free language of differential forms [43] instead of the vector field presentation.Thus, by applying exterior calculus (see, e.g., [20,29,54,61]), we concentrate on the discrete exterior calculus (DEC) of differential forms.The approach follows the groundwork presented by Marsden and his group [23,13,14] and pioneered in electromagnetics simulations by Bossavit and Kettunen [8].Within this framework, the differential operators can be presented exactly at the discrete level.The physical characterization of the discretization is presented by the discrete Hodge operators defining the connection between the primal and dual forms.The discretization is independent of the coordinate systems, and the only source of the discretization error is in the Hodge operator.With the orthogonality of the primal and dual mesh elements, we ensure the diagonal construction of Hodge operators that provides a significant savings in computing time.
The DEC-based discretization in polyhedral meshes was recently studied in [22,45], but now we concentrate on partly structured nonuniform grids, based on the structure of elemental metals and intermetallic compounds.Such natural crystals form stable structures that minimize the potential energy.The structures are constructed from particles packed as close as possible to each other.The feature provides an inspiring starting point for generating spatially isotropic grids that maintain the uniform physical properties in all directions.This is important in electromagnetic simulations in which the waves can propagate in any direction.
Time discretization has a significant role in the efficiency of the time-dependent simulations.Since we use nonuniform grids, we need, instead of conventional staggered leapfrog time discretization [62], a more flexible scheme that adjusts to the grid with varying time step sizes.Lew et al. [33] presented for elastodynamics asynchronous time stepping that was applied to electromagnetic problems by Stern et al. [50].Since the approach needs information from several earlier time steps, we proceed another way by introducing a nonuniform leapfrog method, in which only the quantities at the previous time step are needed to evaluate the variables.
The novelty of the present work is that we present, in section 2, nature-inspired nonuniform grids that have not been used previously in the DEC context.Further, in section 3, we provide an optimized Hodge operator for propagation of plane waves and asynchronous leapfrog-style time discretization that is exact for time-harmonic simulations.In section 4, we validate this state-of-the-art approach with numerical experiments.Concluding remarks are presented in section 5.
2. Discretization mesh.We discretize the three-dimensional spatial domain with a pair of (primal) mesh and dual mesh elements.The primal mesh elements are p-cells with dimension p = 0, 1, 2, 3 as follows: A 0-cell (node) is a vertex v ∈ R 3 .A 1-cell (edge) is defined as a line segment between two 0-cells, and it is oriented from the first to the second node of the list.A 2-cell (face) is a convex planar polygon surrounded by a finite set of edges.An oriented 3-cell (body) is a convex polyhedron surrounded by a finite set of faces with outward-pointing normal vectors.To each p-cell, we assign a dual (3 − p)-cell.
For constructing the dual mesh, we concentrate on the weighted circumcenters to obtain the orthogonality of the primal and dual mesh elements.In practice, a mesh and its (weighted) circumcentric dual mesh can be constructed by a (weighted) Voronoi diagram from a set of (weighted) primal nodes [58,2,11,46].Each Voronoi cell, centered by a node, corresponds to a dual body element.The Voronoi planes (boundaries between the Voronoi cells) correspond to dual face elements, which are assigned to primal edge elements between two nodes.Each Voronoi edge (the intersection of the Voronoi planes) generates a dual edge and a primal face.The dual nodes and the primal bodies are generated from the crossing points of the Voronoi edges.In general, the primal mesh is not simplicial.Convex polyhedral elements, such as cubes, prisms, pyramids, octahedra etc., are possible as well.All of these elements fit easily into the DEC framework.
There are also methods for minimizing the bias between the primal and dual element positions.For instance, Hodge-optimized triangulation (HOT), tailored for the DEC framework (see [40]), can be applied to unstructured triangulations.The idea is to optimize the position and weight of each primal node with an iterative algorithm by keeping the topological relations of the mesh unchanged.
We concentrate on a partly structured grid design and generate a mesh consisting of structured and unstructured (e.g., curved boundaries) areas.This choice combines the geometric flexibility of unstructured grids and the fast generation of structured grids.Since the electromagnetic waves can propagate in any direction, it is important to maintain uniform physical properties in all directions by using spatial grids that re as isotropic as possible.Next, we recall certain natural crystal structures that can be applied with the DEC framework.
2.1.Cubic crystal systems.Cubic tiling is the simplest and most common way to fill a three-dimensional space.It is also the basis of cubic crystal systems, which are used in crystallography to explain natural crystal structures.There are three main structures in cubic crystal systems: primitive cubic, face-centered cubic (FCC) and body-centered cubic (BCC) [10].The primitive cubic structure leads to a regular grid, which has a regular dual grid, corresponding to the conventional Yee scheme [62].
The FCC structure can be constructed by the Voronoi diagram from the vertex set, where the cubic grid vertices and the center points of each face element are united as shown in Figure 1.The space tiling of the FCC structure is composed of alternating regular octahedra and regular tetrahedra in a ratio of 1:2, as shown in Figure 2(a).Each primal edge is of equal length, and the primal faces are equilateral triangles of equal areas.The circumcentric dual faces are rhombuses (parallelograms with four sides of equal length) with relation √ 2 between the longer and shorter diagonal lengths.The circumcentric dual body elements are Kepler's rhombic dodecahedra consisting of 12 congruent rhombic faces [26].
The BCC structure is a tetrahedron structure, which is one of the Sommerville grids [49].The structure is constructed from the cubic grid by adding a vertex at each body center (see Figure 1).The space tiling of the BCC grid is composed of congruent tetrahedra, where each face is a congruent isosceles triangle (see Figure 2(b)).The bottom edges are longer than the side edges by the relation 2 : √ 3. The dual space tiling consists of truncated octahedra; each has six square faces and eight regular hexagonal faces, where all the edges have the same length.This tessellation is known as the Kelvin structure, which was introduced in [55].According to the Kelvin  conjecture, such a polyhedron was proposed as the best space tiling for minimizing the total surface area of the interfaces between elements of equal volume.The BCC structure (or Sommerville's grid) is the simplest well-centered tetrahedral grid for Yee-like simulation schemes, and it is preferred in several instances (see, e.g., [7,57]).Nevertheless, this construction is not likely to be optimal.The dihedral angle of a regular tetrahedron is arccos 1  3 , and thus, in a hypothetical space tiling configuration, we would need 2π/ arccos 1  3 (approximately 5.1) regular tetrahedra to share the same edge [52].From this, it can be concluded that in densely packed nearly regular tetrahedral tiling five or six tetrahedra should share the same edge.Therefore, the dual elements should have either pentagonal or hexagonal faces.This is not the case with either Kepler's or Kelvin's structure.Another class of tetrahedral space tilings has this property, and they are known as tetrahedrally close-packed (TCP) structures.

Tetrahedrally close-packed structures.
A mesh with sharp dihedral angles is preferred when a Yee-like scheme is applied on tetrahedral meshes [7].To achieve grids with such angles, we recall the concept of TCP structures described by Frank and Kasper [18].In all of these structures, there are four combinatorial types of Voronoi cells, which all have only pentagonal and hexagonal faces, with no adjacent hexagons [51].In the literature, the most common TCP structures are A15, C15, and Z lattices.The other known TCP structures can be viewed as combinations of the three basic structures [30,15,48].The maximum dihedral angles of the A15, C15, and Z elements are relatively sharp compared to the BCC structure [17].However, it is not obvious that the TCP structures are better for Yee-like schemes.For example, the edge lengths differ less in the BCC tetrahedra than in the TCP structures [51].The A15, C15, and Z structures are illustrated in Figure 3, and the vertices of each TCP structure are listed in Table 1.
The A15 structure was first discovered in a molecule structure by Hartmann Ebert, and Bretschneider [21].It can be constructed by starting with a BCC lattice
All the grids considered in this paper, except the Z structure, are symmetric in positive and negative x-, y-and z-directions.The Z grid, defined as A 4 B 3 by Frank and Kasper [18], has different properties in all three coordinate directions.Still, the grid is symmetric on the x-y plane in 60-degree increments.The fundamental structure of the Z grid is repeated in h increments in y-and z-directions and in √ 3h increments in the x-direction (see Figure 3).The listed x-coordinates of Table 1 should then be multiplied by √ 3 to obtain the correct tiling.The replicable dual structure consists of six 12-hedra, four 14-hedra and four 15-hedra [51].
3. Computational model.The theory of electromagnetism can be expressed, in a natural way, with the exterior calculus of differential forms (see, e.g., [13,23,50]).The vector fields in model (1.1)-(1.3)are replaced by discrete counterparts such that E = (E 1 , . . ., E n ) T and H = (H 1 , . . ., H m ) T are discrete 1-forms, and T are discrete 2-forms.The discrete 1-form values E i and H j are defined as the integrals of the corresponding vector field over the primal and dual edge elements E i and E * j , respectively.Similarly, the 2-form values are defined as the integrals over the corresponding primal face F j and dual face F * i such that The m × n incidence matrix d 1 represents the neighboring relations and relative orientations of the primal edges and faces.The nonzero entry (d 1 ) j,i = ±1 means that the edge E i is included in the boundary of the face F j .The sign of the entry depends on the relative orientation defined by the counterclockwise circulation.Since , where d 1 and its transpose d T 1 represent the discrete exterior derivatives.
The constitutive relations and the source functions are applied by the concept of the discrete Hodge star operator , which is a mapping from the discrete 1-form to a discrete 2-form.Basically, we determine n × n matrices , σ and m × m matrices µ and σ * to satisfy the equations The spatially discrete Maxwell's system is then presented as where E and H are the vectors of degrees of freedom, associated with the primal and dual edges, respectively.
The construction of the discrete Hodge operator relies on the grid quality.Since we assume the orthogonality of the primal and dual grids, the isotropic material parameters lead to diagonal Hodge operators, which can be presented as diagonal matrices.We apply the diagonal Hodge operators for anisotropic material parameters, but the duality restriction is slightly stronger, as discussed in [44,Section 3.3].The local material parameters are obtained from the tensor fields by where the unit direction vector of the edge E is n E and the convex hull of the face F and the edge E is presented as F E, the volume of which is |F E|.Conventional Yee's discrete Hodge operators [6,14,16,3,40] are based on the assumption of linear vector fields, and the components are expressed as where |F| represents the area of the face F, and |E| is the length of the edge E.
3.1.Harmonic Hodge operator.The only source of the discretization error is in the Hodge operator, the accuracy of which can be enhanced by making assumptions about the nature of the wave and minimizing the errors in (3.3).We concentrate on the special case of time-harmonic plane waves and present the components of the direction-dependent discrete forms E d and D d such that where d is the propagation direction, E d 0 is a complex-valued initial vector, i is the imaginary unit, ω is the angular frequency and p is the spatial position.We seek that minimizes the squared norm of the error First, we consider a coordinate-oriented cubic grid of edge length h, where the primal edge E i is an origin-centered line segment on the z-axis and the corresponding dual face F * i is an origin-centered square on the x-y plane.By assuming a plane wave propagating in homogeneous material in the x-direction with wavelength λ, the direction-dependent exact Hodge operator d is written as (3.9) , where (E d 0 ) z is the z-component of the initial vector and Y ee i,i is the value obtained by Yee's Hodge operator.If the edge length is small compared to the wavelength, the two Hodge operators are close to equal.The larger the elements, the more Yee's Hodge operator overestimates the material terms.Since the systematic error is repeated throughout the domain, Yee's Hodge operator produces a wavelength that is too short.
To obtain a harmonic Hodge operator for more general usage, we integrate the terms of (3.8) in all propagation directions and operate in spherical coordinates with d = (cos(θ), sin(θ) cos(φ), sin(θ) sin(φ)) T .We consider a geometry that is symmetric with the azimuth angle φ and simplify the integration by setting φ = 0. We assume that a primal edge E i , centered at the origin, is pointing in the x-direction, and that the corresponding origin-centered circular dual face F * i , of radius r, is orthogonal to the edge on the y-z plane (see Figure 4).To expand the exponential function in a Taylor series, the edge and the face element are assumed to be small compared to the wavelength.
We concentrate on the x-component of the initial vector (E d 0 ) x and set an auxiliary variable α := iω √ i µ i sin θ.We also change the variables such that z = r sin u, dz = r cos u du, and √ r 2 − z 2 = r cos u to obtain the Wallis integral (see, e.g., [47]), and we get the power series presentation Respectively, we set an auxiliary variable β := iω √ i µ i cos θ and formulate E d i by onedimensional integration over the edge element to get the power series presentation By minimizing the error R d i in the squared norm with i,i = harm i,i , we get the harmonic Hodge operator where κ i is the curvature correction for the Hodge term harm i,i , κ F = ω 2 i µ i r 2 , and The harmonic Hodge operators for magnetic permeability, electric conductivity, and magnetic conductivity can be presented as, respectively, (3.12) where the curvature correction κ * j is computed similarly to (3.11), but between the elements F j and E * j instead of F * i and E i .From the computational point of view, the dual faces are not circular but polygonal.That is why we compute an approximation of r 2 , involved in (3.11), for arbitrarily shaped polygonal face elements.Essentially, the maximum internal sphere of radius r min and the minimal external sphere of radius r max define the lower and upper limits for the radius of the element (see Figure 5).For regular polygons, r min corresponds to the distance between the face center and the midpoint of the edge.Respectively, r max corresponds to the distance between the face center and the vertex.By using r min r max  the approximation r 2 ≈ 1 3 (2(r min ) 2 + (r max ) 2 ) for deriving κ i and κ * j , the first three terms of the Taylor polynomials in (3.11) are exactly those obtained by the polynomial integration.In general, for an irregular polygon with n edges, we approximate r 2 as a mean value (2(r min j where r min j and r max j are the internal and external radii associated with the jth edge and vertex (see Figure 6).We approximate the edge length as the average of the n edge lengths.Obviously, the number of operations needed to compute r 2 for an irregular element is n times that for a regular element.

Time discretization.
Typically, staggered leapfrog time discretization is used in the Yee-like schemes due to its simplicity and energy-preserving properties (see, e.g., [62,13,53,60,36]).The method is only of second-order accuracy, but it can be improved by using additional information about the time evolution of electromagnetic waves.In the FDTD context, there exists an exact leapfrog method for time-harmonic problems [35].In what follows, we present a modification, generalized to the DEC framework.We also proceed further and propose an entirely new nonuniform version of the method.
First, we set t 0 as the initial time and divide the time period T into time steps, each of which has length ∆t.At the kth time step, the time-discretized form of the variables E and H are defined as where t k = t 0 + k∆t.Since we consider time-harmonic problems, we set E(t) = Re( Êe −iωt ) and H(t) = Re( Ĥe −iωt ), where Ê and Ĥ are complex-valued spatial variables.This presentation leads to the following exact time discretization formulas: Time stepping is obtained by starting with the initial conditions E 0 and H 0 and computing E at time t k + ∆t 2 and H at time t k+1 , k = 0, 1, 2, . . ., (see Figure 7), by Basically, the formulas are obtained from conventional Yee's leapfrog method by scaling the discrete Hodge operators and µ by ϕ sin ϕ and σ and σ * by 1 cos ϕ , where ϕ = ω∆t 2 .Thus, the harmonic leapfrog method has the same convergence and energy conservation properties as Yee's leapfrog method.
3.2.1.Stability.The leapfrog method described above is conditionally stable.That is, the stability of the time discretization depends on the length of time step ∆t, which, for its part, must be adjusted to the grid size.Since we use, in general, unstructured grids, we optimize the time step size for each element and validate the stability condition.
We start by formulating a stability criterion that prevents excessively large changes in the differences between the successively computed values.We assume a nonnegative real value E max , limiting the values of E i such that |E i | ≤ E max i for all i (see Figure 8), implying that That is, we can choose a constant C that defines, by the formula how large the allowed changes can be.If C = 4, then E i is allowed to vary between −E max i and E max i on the consecutive time steps.Thus, we can set an upper limit C ≤ 4.
However, by assuming that σ = σ * = 0 and applying (3.15)-(3.16)with the time step size ∆t Ei for E i , we can see that By combining (3.17) and (3.18), we can write the stability criterion as   In nonuniform time stepping, ∆t is divided into smaller steps only where this is needed.In this 1-dimensional illustration, the circles with a number inside represent time instances of E i and H j .The arrows illustrate the chronological order of the computation.
Since (2/ω) sin(ω∆t Ei /2) < ∆t Ei , we can write a sufficient condition to fulfill the stability criterion (3.17) for E i with the time step size ∆t Ei , and for H j with the time step size ∆t Hj , such that, The choices of E max i and H max j determine the distribution for the nonuniform time step sizes.For the numerical experiments presented in section 4, we use the limit values because we have discovered that this choice enables an efficient simulation.

Nonuniform leapfrog method.
The nonuniform leapfrog method is initialized by computing the maximal local time step sizes ∆t Ei and ∆t Hj by (3.20).After that, the global time step ∆t, which is the smallest uniform iteration block to be repeated, is selected between the minimal and maximal values of the local time step sizes ∆t Ei and ∆t Hj .
The global time step is divided for each element E i and H j by integer numbers s Ei and s Hj , respectively (see Figure 9).For deriving (3.20), we assumed that the values H j , neighbors of the values E i , are discretized by using the same time step size ∆t Ei .For nonuniform time stepping, we loosen the condition by requiring that the neighboring values be discretized at the same level of the time step size.That is, if the primal edge E i corresponding to E i lies on the boundary of the primal face F j corresponding to H j , we set These inequalities offer a small enough time step size for each element, but the energy conservation is not guaranteed, since the time stepping is asynchronous.We  4. Numerical experiments.In this section, we report the numerical experiments for the accuracy of the Hodge operator using different grids.We start the experiments by considering the stability of the nonuniform leapfrog method.After that, we concentrate on the simulated wavelength in different wave propagating directions.The last numerical example shows how the wavelength error is correlated to the overall efficiency of the method.All simulations are run with the grids presented in section 2, and the properties of each grid are analyzed.Yee's Hodge and harmonic Hodge approximations of section 3 are also compared.The algorithms are implemented in the C++ programming language.
4.1.Stability of the nonuniform leapfrog method.In the first example, we focus on the stability, computational efficiency, and energy conservation of nonuniform time stepping.The simulation domain Ω is a cube of edge length 2 that is discretized using the randomly generated mesh (see Figure 10(a)).The boundary of the mesh is constructed with squares of edge length 0.1.The interior nodes are generated by a uniform density function resulting in mesh density of one node per volume 10 −3 .The minimal distance between nodes is controlled by regenerating each node that is less than 0.01 from another node.The primal and dual meshes are constructed with the Delaunay triangulation method, where the three-dimensional Voronoi diagram is applied [12,38,41,2,11].The range of the element sizes is wide, as shown in Table 2.We solve (1.1)-( 1.3), where the material parameters are set as = µ = 1 and σ = σ * = 0.At t 0 = 0, the initial terms E 0 = (0, cos (2πx) , sin (2πx)) T and H 0 = (0, − sin (2πx) , cos (2πx)) T satisfy a circularly polarized plane wave with wavelength λ = 1 propagating along the x-axis, through the y-z plane.The wave is trapped inside the domain by setting on the boundary of the domain the Dirichlet boundary condition, E×n = 0 and H•n = 0, where n is the outward-pointing unit normal vector.For the time interval [0, 1000T ], where T = 1 is the period of the electromagnetic wave, we apply uniform and nonuniform time stepping schemes with stability constant C = 4.The minimum and maximum values of the time step factors s Ei and s Hj in (3.23)-(3.24)are 1 and 81, respectively, which means there is a remarkable difference between the smallest and largest time step sizes.In the uniform time stepping scheme, we use the smallest time step involved in the nonuniform scheme.We performed the simulations on a single Intel Xeon E5-2670 processor at 2.60 GHz.The CPU time used for one time period T with the uniform leapfrog method was 17.8 seconds, while the time consumption with the nonuniform leapfrog method was 1.05 seconds per time period T .Thus, the nonuniform time stepping caused significant improvement in the simulation efficiency.
The energy P := 1 2 E T E + H T µH was computed after each time period and is presented in Figure 11.The numerical experiments show that the energy conservation is exact with the uniform leapfrog method.When the nonuniform leapfrog method is used, the energy varies slightly with time.However, the values of the energy seem to be bound between upper and lower limits.The standard deviation of the energy is 0.0040, which means the relative standard deviation compared to the mean value is 0.028%.Thus, we can deduce that the energy is conserved in the long term.
The results indicate that the system is stable with uniform and nonuniform time stepping schemes with stability constant C = 4.With another simulation, we found that neither time stepping scheme produced a stable system with constant C = 4.5.Thus, the stability criterion nominates a very tight bound for the time step size in the current random mesh discretization.However, we have found that with certain spatial discretizations the constant C = 4 is not sufficient for the convergence.Thus, we prefer using term C = 2 to obtain a reliable criterion, and this constant is applied in the simulations.

Simulated wavelength.
In the second experiment, we consider the simulated wavelength in different propagation directions.The computational domain  is a sphere of radius 5. Inside the domain, we use the same parameters as in section 4.1.The electromagnetic wave propagates with the angular frequency ω = 2π at speed c = 1; that is, the theoretical wavelength is λ = 1.To let the electromagnetic waves travel out from the computational domain without reflections, we surround the domain with an absorbing boundary layer [24] of thickness 2λ (see Figure 12).The absorbing layer is set such that = µ = 1, and the electric and magnetic conductivities, σ and σ * , satisfy the condition σ = σ * µ = 0.4Υ, where Υ is the outward distance from the surface of the inner sphere.The simulation is started with the initial values E 0 = 0 and H 0 = 0, and the circularly polarized incident waves E inc = (0, cos (2π(x − t)) , sin (2π(x − t))) T and T are generated by a source function defined on the absorbing layer.
For spatial discretization of the domain and the absorbing layer, we use six different grids, presented in section 2. These are cubic, FCC, BCC, A15, C15, and Z grids.The lengths of the edge elements in the cubic grid are h = λ 10 , whereas the other grids are scaled to keep the computing time fixed (see Table 3).To validate the spatial isotropy, the simulations are run in 10 different wave propagation directions, d, labeled from A to J (see Figure 13).All the grids, except the Z grid, are symmetric in positive and negative x-, y-, and z-directions.With those grids, the selection of directions is comprehensive.The directions from A to J cover 1/48 of the full unit   sphere area, and due to overlapping directions, the symmetry directions of the 10 cases cover 218 directions in the full unit sphere.The Z grid is tested in six different coordinate axis permutations to cover all these directions.The angle between the adjacent directions is about 15 degrees.
The grids are uniform throughout the domain, implying that the time stepping scheme (3.23)-(3.24) is applied with s Ei = s Hj = 1.Thus, the time stepping strategy is equal to the uniform leapfrog method (3.15)- (3.16).For each wave propagation direction, the values of the simulated wavelength are observed after the simulation is run through 500 time periods on a line of length 8λ, which is parallel to the wave propagation direction and centered at the domain center (see Figure 12).The phase error of the simulated wave is computed at several positions on the line.The position and the phase error are plotted on an x-y grid, and a simple linear regression curve is fitted on the result.The slope defines the wavelength error.Yee's Hodge and the harmonic Hodge operators are considered, and the errors of the directional wavelengths are illustrated in Figure 14.The relative L 2 norm of the difference of the solutions at the end of the last two time periods was less than 10 −6 ; in all cases, this implies a well-converged result.
The wavelength error statistics are reported in Table 4.In the case of Yee's Hodge operator (see Figure 14(a)), the largest error, 1.63%, occurs with the cubic grid in orthogonal case A. This error seems to be close to the factor between the exact Hodge operator d and Yee's Hodge Y ee in (3.9).Naturally, the computed wavelength varies with the propagation direction.On average, the most accurate results are achieved with the FCC grid.In that case, the average error is 0.82%, while for the other grids, it is between 0.96%and 1.00%.However, the tetrahedral grids have the smallest directional dependencies for the wavelength.The cubic grid has the largest directional dependency with a standard deviation of 0.37%.On the whole, the same characteristics are observed with the harmonic Hodge operators (see Figure 14(b)).However, the average wavelength is moved very close to the exact solution (error ≤ 0.04%) when the harmonic Hodge operator is considered.It seems that the C15 grid is the most isotropic; with this grid the standard deviation of the error is 0.03%.The current consideration shows how the wavelength varies by the propagation direction in staggered grids.Different grids have unique properties, and the wavelength error in each direction depends on the element size and the selection of the grid.If the harmonic Hodge operator is the systematic error on the wavelength is very small.Next, we will study how these results transfer to the simulation of a scattering problem.

4.3.
Scattering by a sphere.We consider the accuracy and efficiency of the method in a spherical computational domain of radius 2.7, in which the spherical scatterer of radius 2.5 is centered.Inside the scatterer, the material parameters are = 2.5599, σ = 0.064π, µ = 1, and σ * = 0.The material parameters = µ = 1 and σ = σ * = 0 are applied in the rest of the computational domain.Accordingly, the wavelength is λ = 0.625 inside the scatterer and λ = 1 in the computational domain outside the scatterer.
The electric and magnetic fields can be considered sums of the incident waves E inc and H inc and the scattered waves E sca and H sca , such that E = E inc + E sca and H = H inc + H sca .Thus, we concentrate on the accuracy and efficiency of the scheme by solving (1.1)-(1.3)for the scattered wave fields E sca and H sca .The incident wave components are E inc = (0, cos (2π(x − t)) , 0) T and H inc = (0, 0, cos (2π(x − t))) T , presenting a fully polarized plane wave of angular frequency ω = 2π and time period T = 1, and propagating in the direction of the positive x-axis.
Essentially, we use the same grid types and Hodge operators as in the experiments presented in section 4.2.The interior of the scatterer is discretized by using four different discretization levels for each of the six grid types.The mesh resolution for the simple cubic mesh is defined such that the simulations are run with the edge element lengths 1 11.2 , 1 16 , 1 22.4 , and 1 32 , corresponding to 7, 10, 14, and 20 elements per wavelength inside the scatterer.The elements of the other grids are scaled by the factors presented in Table 3.
The boundary of the scatterer is covered by triangles (see Figure 15(a)) and the space between the boundary and the interior grid is constructed with the Voronoi tessellation.The boundary triangulation is optimized by the HOT optimization method, discussed in section 2, to improve the element quality.Further, the boundary surface is stretched in the radial direction to the domain surface, so that a layer of  thickness 1.7, filling the domain and absorbing layer, is generated outside the scatterer.The elements outside the scatterer are scaled to match each discretization level (see Figure 15(a)).The 1.5 thick perfectly matched layer (PML) [4] is set outside the domain using the conductivities σ and σ * for the outgoing wave by the relation σ = σ * µ = 0.5Υ, where Υ is the outward distance from the surface of the inner sphere.An asynchronous time stepping scheme with harmonic time stepping formulas was applied in these simulations to obtain an optimal time step size for each part of the mesh, as presented in section 3.2.2.
For comparison, we also perform the same simulations with the classical Yee scheme, where the simulation domain is discretized by uniform cubic elements.The spherical scatterer of radius 2.5 is centered in a cubic domain, of edge length 5.4, corresponding to the smallest cube that can contain the spherical domain of radius 2.7 (see Figure 15(b)).To maintain the same time consumption as in the other scattering simulations, the edge element lengths 1  8 , 1 11.2 , 1 16 , and 1  22.4 are applied in the whole simulation domain.The domain is surrounded by the PML of thickness 1.5 with the same formula as in the case of the other grids.Yee's Hodge operator and Yee's leapfrog time discretization scheme are applied.
In the beginning of the simulation, the variables were initialized with the values E 0 = 0 and H 0 = 0.After 500 time periods were simulated, the relative L 2 norm of the difference between the consecutive periods decreased to less than 4 • 10 −4 , in each case, implying a relatively well-converged solution.To elaborate the accuracy of the method more precisely, the simulated fields are observed on the boundary between the domain and the PML.Then, the near-field solution is transferred to the far-field solution by the technique applied to the FDTD method in [56].The farfield scattering data is applied to produce the Mueller matrix [5], which includes the scattering intensities and polarization in all scattering directions.Figure 16 illustrates the relative error, measured by the L 2 norm, of Mueller matrices compared to the exact Mie solution [9] computed by the Mie scattering code presented in [37].
The simulation time increases as the element size becomes smaller.Naturally, a smaller element size implies a finer discretization and a smaller error.For each grid type, the relative error is smaller in the case of the harmonic Hodge operator compared to Yee's Hodge operator.The simple cubic mesh with the harmonic Hodge operator seems to produce about equally efficient simulations as the FCC grid with Yee's Hodge operator.The grid types affect the accuracy in a similar way as expected by the results presented in section 4.2.With Yee's Hodge operator, the FCC grid produces the most accurate results.By using the FCC grid, a particular accuracy, expressed as relative error, is achieved with considerably lower CPU time than by using any other grid with Yee's Hodge operator (see Figure 17).For instance, if we set the target relative error to be 5-10%the CPU time required to solve the problem with the FCC grid is around one fifth of that required with the simple cubic mesh and less than one tenth compared to the Yee approach.As a matter of fact, all of the bars presenting the Yee scheme are truncated from the top in Figure 17.That is because the scheme needs more than 7500 seconds to provide any of the considered accuracies.This is also the case with the schemes related to the other truncated bars.
As time-harmonic problems are considered, even higher CPU time savings are gained by utilizing the harmonic Hodge operator.In that case, the simulation with the simple cubic mesh demands the most CPU time, while the most accurate results are obtained with the C15 and Z grids.The simulations with the C15 grid reach a prescribed accuracy five times faster than the simulations with the simple cubic mesh.Further, the simulations based on the harmonic Hodge operator can be even faster with the Z grid than with the C15 grid if high accuracy is required.Thus, by selecting an appropriate grid and a sophisticated Hodge operator, we managed to decrease the time consumption to a small percentage compared to the conventional methods based on cubical grids and Yee's Hodge operator.

5.
Conclusions.We considered time-dependent electromagnetics in the framework of spatial discretizations based on quasi-uniform three-dimensional polyhedral grids, inspired by the structure of natural crystals and by discrete exterior calculus.By applying the conventional leapfrog time integration, we obtain a general purpose scheme for transient simulations in complex media.Since there is a wide range of application areas, such as remote sensing and terahertz spectroscopy, that benefit from time evolution of time-harmonic waves or pulses, we further improved the accuracy of the approach for time-harmonic waves.We optimized the Hodge operator for decreasing the error of the simulated wavelength and presented an asynchronous leapfrog-style time discretization that provides enhanced efficiency with non-structured grids.We showed that by selecting an appropriate grid and a sophisticated Hodge operator, the simulation can be run using only a small percentage of the computing time needed with the classical Yee scheme.

Figure 1 .
Figure 1.The vertices of cubic, FCC, and BCC grids are represented by the • symbol on even layers and by the + symbol on odd layers.

Figure 2 .
Figure 2. The FCC and BCC structures and the corresponding regular grids.The vertices of the primal and dual elements are represented by small spheres.

Figure 3 .
Figure 3.The vertex positions and primal and dual structures of the TCP structures A15, C15, and Z.

Figure 5 .
Figure 5.The internal and external spheres define the lower and upper limits of the radius r of the element (dotted sphere).

4 Figure 6 .
Figure 6.For an irregular quadrilateral, r is replaced by the mean value computed from the distances r min j and r max j , j = 1, . . ., 4.

Figure 9 .
Figure 9.In nonuniform time stepping, ∆t is divided into smaller steps only where this is needed.In this 1-dimensional illustration, the circles with a number inside represent time instances of E i and H j .The arrows illustrate the chronological order of the computation.
(a) Random mesh.(b)The initialized wave is trapped inside the domain.

Figure 10 .
Figure 10.The energy conservation is simulated in a cube of edge length 2.

Figure 11 .Figure 12 .
Figure 11.Uniform time stepping conserves the energy exactly.With nonuniform time stepping, the energy remains almost constant during the long-term simulation.

Figure 13 .
Figure 13.Directions of simulation cases illustrated on a cross section of a unit sphere.

Figure 14 .
Figure 14.The relative error of the simulated wavelength.Each grid is tested with 10 propagation directions, labeled from A to J, and the non-symmetric Z grid is additionally rotated in six different orientations to cover all the directions.
(a) The spherical domain with a cubic grid inside the scatterer.(b) The domain discretized by cubes for the Yee scheme simulation.

Figure 15 .Figure 16 .
Figure 15.Examples of discretizations for simulating scattering by a sphere.

Figure 17 .
Figure 17.CPU times needed to reach target accuracies.

Table 2
The element sizes of the random Voronoi cube.Ei and s Hj such that they can be written as a power of three, i.e., s = 3 u ∈ {1, 3, 9, 27, ...}.With this selection, each neighboring element, which has finer time stepping than the current element, also has a midpoint instance for the current value update (see Figure9).This property is the basis for the leapfrog methods, and thus, we call this scheme a nonuniform leapfrog method.The time stepping of the nonuniform leapfrog method is carried out in chronological order.The value updates are strictly iterative, meaning that a new value E new i or H new j can be computed from the latest instance of values.Thus, no extra storage is required during the computation.The following two equations illustrate the value updates in the nonuniform leapfrog method, where the new values replace the previous values immediately.Similarly to (3.15)-(3.16),we can write

Table 3
Primal and dual (*) element sizes of scaled grids.The grids are scaled by constants 1.0, 2.12, 1.84, 2.92, 4.20, and 2.92.The corresponding computing time for one time period with h = 0.1 is reported in seconds.

Table 4
Statistics of the wavelength error: average, standard deviation, and minimum and maximum of the absolute values.