High-quality discretizations for microwave simulations

We apply high-quality discretizations to simulate electromagnetic microwaves. Instead of the vector field presentations, we focus on differential forms and discretize the model in the spatial domain using the discrete exterior calculus. At the discrete level, both the Hodge operators and the time discretization are optimized for time-harmonic simulations. Non-uniform spatial and temporal discretization are applied in problems in which the wavelength is highly-variable and geometry contains sub-wavelength structures.


I. INTRODUCTION
Microwaves are high-frequency electromagnetic waves that can propagate in any spatial direction. Thus, it is important to maintain uniform physical properties in all directions, which is not the case when the cubic grid is applied in the finitedifference time-domain (FDTD) method, originally proposed by Yee [1]. The quality of the discrete model, based on the structure of the computational grid, has a significant impact on the efficiency of the method. We use non-uniform grids and associate the degrees of freedom of the electric and magnetic fields with primal and dual grids, respectively. The physical characterization of the discretization is presented by the discrete Hodge operators defining the connection between the primal and dual forms [2]. The discretization is independent of the coordinate systems, and by the orthogonality (with respect to metric) of the primal and dual mesh elements, we ensure the diagonal construction of Hodge operators providing significant savings in computing time [3].

II. ELECTROMAGNETICS IN THE TIME DOMAIN
The work needed to move a particle with electric charge q, along an oriented curve C, in an electric field E = (E 1 , E 2 , E 3 ) T is where dl is the differential unit tangent vector of the curve, dx i , i = 1, 2, 3, are differentials and E = E 1 dx 1 + E 2 dx 2 + E 3 dx 3 is a differential 1-form. Respectively, we can consider a magnetic charge and obtain a differential 1-form H = H 1 dx 1 + H 2 dx 2 + H 3 dx 3 (instead of the magnetic field H).
The current flowing through a surface S where dA = n · dS is the differential vector element of surface area A normal to surface S, ∧ is the exterior product (wedge product), and electric current density Respectively, we can define the magnetic current density J * . The electric charge in volume V is where is the electric charge density and Q = dx 1 ∧ dx 2 ∧ dx 3 is a differential 3-form presenting the electric charge. Respectively, Q * = * dx 1 ∧ dx 2 ∧ dx 3 the magnetic charge. The generalized Stokes theorem, applied to a differential form A over the boundary of an oriented manifold Ω, is where d is the exterior derivative. The Hodge star , defining the coordinate system and metric properties, maps a differential k-form to a differential n − k form (n = 3 is the spatial dimension). The constitutive relations are written as D = εE and B = µH, where D and B are differential 2-forms and ε is the electric permittivity and µ is the magnetic permeability. The electromagnetic waves are presented by where the exterior derivative d on 1-forms corresponds to the curl operator on vector fields whereas on 2-forms it corresponds to the div operator on vector fields. By taking the exterior derivative of Equations (5) and (6), we get the continuity of the charges, and the information given by Equations (7) can be covered by the initial conditions. III. DISCRETIZATION The computational domain Ω is discretized by a mesh, in which we have vertices, edges, surfaces, and volumes. In particular, we are interested in the edge elements E j , j = 1, . . . , N E and face elements F i , i = 1, . . . , N F , where N E and N F are the number of the edge and face elements, respectively. We also construct a dual mesh with edge elements E * i , which are duals to the primal face elements, and face elements F * j , which are duals to the primal edge elements.

A. Discrete Exterior Calculus
The spatial discretization is applied using the discrete exterior calculus (DEC) [4]. The discrete 1-form E and the discrete 2-forms B and J * are presented as column vectors of components The discrete exterior derivative d (boundary operator) is a N F × N E matrix presenting the incidence number of each the value at the ith row and the jth column of the incidence matrix is ±1, depending on the relative orientation between the face and the edge element. The spatially discretized formulation of Equation (6) is On the dual mesh, we have discrete dual forms. The discrete dual 1-form H and the discrete dual 2-forms D and J are presented as column vectors of components The elements of the dual mesh have a similar orientation to the corresponding primal elements if the incidence number of each E j to each face F i equals to the incidence number of each E * i to each face F * j . Then, the dual discrete exterior derivative is d T and the spatially discretized formulation of Equation (5) is presented as The discrete Hodge stars are needed to present the constitutive relations D = εE and B = µH. The ε-related discrete Hodge ε is a N E × N E diagonal matrix, presenting the linear map from the discrete primal 1-form E to the discrete dual 2-form D. The diagonal terms are where the curvature correction κ j is set to compensate for the systematic wave speed error of harmonic wave. As derived in [3], where κ E = ω 2 µ|E j | 2 and κ F = ω 2 µ(2r 2 + R 2 )/3, and r and R are the inner and outer radii of dual face F * j . The µ-related discrete Hodge operator µ is a N F × N F diagonal matrix mapping the discrete dual 1-form H to the discrete primal 2-form B presented as an N F × N F diagonal matrix where the curvature correction κ * i is derived in a similar manner with respect to the elements F i and E * i . Applying DEC on a regular grid without curvature correction results in FDTD discretization as a special case. The curvature correction and high-quality meshing can increase the method efficiency by orders of magnitude [3].

B. Time Discretization
The time evolution begins with the initial conditions E 0 = E(− ∆t 2 ) and H 0 = H(0), after which it proceeds, as a staggered procedure, by computing E k = E(k∆t − ∆t 2 ) and H k = H(k∆t), where k = 0, 1, 2, . . . and ∆t is the length of the time step. Since we are considering time-harmonic simulations, we optimize the time discretization for a single wave frequency. That is, instead of the conventional leapfrog time-discretization, we apply the formulas where s is the time step size factor [3]. For uniform grids, we can choose s = 1 and apply the CFL condition to get the upper bound for ∆t.
Since we use non-uniform space discretization, and the time step size needs to be adjusted to the wavelength, we allow smaller time steps for the finer parts of the grid. The local time steps are obtained from the localized CFL-type condition, which take only a small neighborhood of elements into account [3]. The time step is adjusted by the integer-valued factors s Ej and s Hi , such that, for updating the values of E k j , we use s = s Ej and for updating the values of H k i , we use s = s Hi . To guarantee the energy conservation, we synchronize the time stepping by setting s Ej , s Hi ∈ {1, 3, 9, 27, . . . }. Using the time stepping approach, we obtain an explicit and causal simulation method, which can be efficiently applied to solve time-harmonic problems [3].

IV. NUMERICAL EXPERIMENTS
The DEC discretization can be applied for arbitrary polyhedral meshing with orthogonal duality. The dual mesh is obtained as a Voronoi diagram generated for a given vertex set and given metric induced by the material. The primal mesh is determined by the topology. The numerical experiments of this section exploit the body-centered cubic (BCC) structure [5], [6] defining the primal grid with congruent tetrahedra and the dual grid with truncated octahedra [7]. A uniform grid is rarely optimal discretization for numerical simulations. Thus, we apply structured grid refining, where vertex structures are applied in different scales. In such a manner, one can efficiently obtain non-uniform grids, where element properties are guaranteed. We apply grid refinement due to the highly variable material parameters and for the discretization of subwavelength structures.

A. Highly variable material parameters
We consider a microwave oven, of which the dimensions are 28.0 × 28.0 × 21.0 centimeters representing depth x, width y, and height z, respectively. Placed in the middle of the wall facing in the positive y-direction, a square-shaped magnetron with an edge length of 12.2 cm produces a harmonic wave of frequency f = 2.45 GHz, implying that the wavelength is λ = 12.2 cm. Inside the microwave oven, centered at height of 3 cm above the bottom of the computational domain, is a bowl, having the shape of a conical frustum. The bowl is filled with liquid with a depth of 6 cm; the bottom radius is 3 cm and the top radius is 6 cm. We approximate the liquid by permittivity ε = 70ε 0 , where ε 0 is the air (or vacuum) permittivity. The absorption term of the liquid is determined by the right-hand side source term J = 0.58ωε 0 E. The structured grid refining was applied involving at least 11,000 unknowns per λ 3 and a total of 2,005,708 unknowns in the domain (see Fig. 1). The perfectly reflecting walls are defined by the Dirichlet boundary condition. The magnetron is modeled by the absorbing boundary condition with source terms E 0 = (cos ωt, 0, sin ωt) T , H 0 = (sin ωt, 0, − cos ωt) T , where ω = 2πf is the angular frequency.
The numerical experiment was initialized by E 0 = 0 and H 0 = 0, and was continued until time t = 501T was reached, where T = 1/f is the wave period. During the first period, the source terms were increased from zero to full terms by the transitions suggested in [8]. For the last 500 iterations, the full source terms were applied. The time step size factor s varied between 1 and 27 and a total of 205,557,312 value updates were executed per iteration. The simulation was carried out in around 4 minutes on 32 Intel Xeon cores at 2.67 GHz. The results, presenting the cross sections of electric and magnetic fields, are illustrated in Fig. 2.
In comparison, an FDTD simulation with ten elements per wavelength means 6000 unknowns per λ 3 . The liquid wavelength is 1.46 cm designating 1940 unknowns per cm 3 . That makes about 31.9 million unknowns in the whole domain. Such a grid has 83.7 elements per vacuum wavelength, and therefore the CFL condition requires 145 time steps per period. In total, the FDTD simulation would take 4.62 billion value updates per period, which is more than 20 times the number of value updates in our numerical simulation.

B. Sub-wavelength structure
In x-direction, the domain is a three-layer structure with thickness 1 truncated by the absorbing boundary condition. The back and front layers are discretized by yand z-periodic mesh, as illustrated in Fig. 3, and the source term for a circularly polarized plane wave of wavelength λ = 1 is set at the back wall. The middle layer is a plane of a periodic sub-wavelength grille structure with circular holes surrounded by perfectly conducting material (E = 0). For the numerical tests, we apply six different grilles, where the periodic pattern is repeated 1, 2, 4, 8, 16, or 40 times per unit length. We apply four discretization levels, which are called 5, 10, 20 and 40. Each discretization has a base grid of the BCC structure, where the maximum edge lengths are 1 5 , 1 10 , 1 20 , and 1 40 , respectively. Structured grid refinement is applied near the plane x = 0 by inserting vertices of the regular structure The simulations were initialized by E 0 = 0 and H 0 = 0 and processed by 10 periods with increasing source terms and by 10 periods with full source terms. The solutions for each discretization level and for the largest grille size are illustrated in Fig. 4. The energy penetration is computed for each grille and each discretization level and the results are illustrated in Fig. 5. The larger the holes are, the more energy penetrates, and for the largest hole size the penetration is approximately 50 %. This is surprising, as holes cover only π 8 ≈ 39 % of the whole surface.
The solutions for each discretization level are close to each other, but upon closer examination, there is up to 25 % relative difference between the results. The conclusion is that the refinement should cover a wider volume. We increase the depth of the refinement such that each regular group of nodes is limited by x ∈ [−kh, kh], where k = 1, 2, 3, 4 is called the refinement thickness. Figure 6 shows that wider refinement improves the accuracy. With the largest refinement width, one has less than 1 % error. The number of unknowns are 1,174,400, 1,466,400, 1,854,400, and 2,146,400 for these cases, respectively, which are less than one tenth compared to the reference discretization.

C. Combined problem
We repeat the microwave simulation by applying a grille with circular holes with diameter 0.22 cm on the wall facing in positive x-direction. The elements size at the grille is 0.11 cm and the discretization is illustrated in Fig. 7. The number of unknowns is 3,729,426 and the time step size factor s vary between 1 and 27. The total number of value updates per period is 751,607,208. The simulation through 501T took about 14 minutes on 32 Intel Xeon cores at 2.67 GHz. As a result, 0.063 % of the power produced by the magnetron is penetrated through the grille. The rest is absorbed as heat in the liquid. The energy loss is so small that the grille does not have a remarkable effect on the progress of electromagnetic energy inside the liquid, as illustrated in Fig. 7. V. CONCLUSIONS We applied polyhedral structures and grid-refining methods with DEC. The grid refinement is useful in modeling tasks with highly variable material parameters and sub-wavelength structures. Further research should consider how the error estimation methods can be applied to find the optimal refinement.