A p-Strategy with a Local Time-Stepping Method in a Discontinuous Galerkin Approach to Solve Electromagnetic Problems

Abstract

We present a local spatial approximation or p-strategy Discontinuous Galerkin method to solve the time-domain Maxwell equations. First, the Discontinuous Galerkin method with a local time-stepping strategy is recalled. Next, in order to increase the efficiency of the method, a local spatial approximation strategy is introduced and studied. While preserving accuracy and by using different spatial approximation orders for each cell, this strategy is very efficient to reduce the computational time and the required memory in numerical simulations using very distorted meshes. Several numerical examples are given to show the interest and the capacity of such method.

1. Introduction

In early works, we have proposed an efficient new Discontinuous Galerkin (DG) method to solve the Maxwell equations in time domain [1]. However, in order to treat industrial problems, we have generally to cope with a very distorted meshes with a large difference between the sizes of the cells. In the DG simulations, the spatial approximation order taken into account is obtained by considering the size of the largest cell and the largest frequency in the spectrum of the excitation source. Consequently, our method suffers from a sampling that is too high for the smallest cells in the mesh and from a time step that is very small to ensure stability. To improve the choice of the time step, we have proposed a local time-stepping strategy in the method [2] which allows reduction of the computational time. Nevertheless, the spatial order of approximation remains the same for all the cells, and for a given source, some cells are too sampled and lead to a local time step on these cells still too small. To avoid this problem, a method with local order on the cells can be used [3]. Some works have been already realized for Maxwell's equation on this subject, with classical upwind DG formulation and ADER strategy [4]. In this paper, we present, for a particular DG method well adapted to the Maxwell equations, a strategy to mix spatial orders of approximation in each cell on the whole mesh. This strategy allows decreasing the order of the small cells and then increasing the CFL condition for these cells and consequently their local time-step. Indeed, it has been shown, for the DG method used, that the CFL condition decreases when the spatial order of approximation increases [5].

In the first part of this paper, we recall the mathematical formulation of the DG schema studied and the local time-stepping strategy used. In the second part, the formulation of a local spatial approximation order or mixed-spatial order strategy is given and introduced in the presented DG scheme. A mathematical study for the stability of the resulting numerical scheme is also shown. Finally, in the third part, a numerical study by using the mixed-spatial order strategy with the local time-stepping strategy is done.

2. Mathematical Formulation

In this section, we recall first the formulation of the DG method taken into account in the study of our mixed-spatial order strategy. Next, we present the local time-stepping strategy used to improve the DG performances.

2.1. Discontinuous Galerkin Method

Let Ω be a bounded open subset of ℝ³ whose boundary is 𝜕Ω, and let 𝑛 denote the unit outward normal to Ω. Let 𝜀(𝑥), 𝜇(𝑥), and 𝜎(𝑥) denote, respectively, the permittivity, the permeability, and the conductivity of the medium.

We consider the problem described by the Maxwell equations.

find (𝐸,𝐻)∶Ω×]0,𝑇[→ℝ³×ℝ³ such that 𝜀𝜕_𝑡𝐸+𝜎𝐸=∇×𝐻,𝜇𝜕_𝑡𝐻=−∇×𝐸,𝑛×𝐸=0on𝜕Ω, (1) where 𝐸 and 𝐻 define the electric and magnetic fields. The boundary condition is not restrictive because it is used both to treat closed problems such as cavities and to bound the Perfectly Matched Layers (PML) domain [6]. Consider a set 𝒯 of hexahedral elements (𝐾_𝑖)_{𝑖=1⋯𝑁} being a partition of Ω. We introduce the following space of approximation: 𝐕_𝑟=𝐯∈𝐿²(Ω)³;∀𝐾∈𝒯,𝐷𝐹^∗_𝐾𝐯|_𝐾∘𝐹_𝐾∈𝑄_𝑟𝐾³, (2) where 𝐾=[0,1]³ is the unit cube or the reference element, for all 𝐾∈𝒯,𝐹_𝐾∶𝐾→𝐾 denotes the trilinear mapping which associates the vertices of each element, 𝑄_𝑟(𝐾) is the space of polynomials of degree at most equal to 𝑟∈ℕ^∗ in each variable on 𝐾, and 𝐷𝐹_𝐾 and 𝐽_𝐾 are, respectively, the Jacobian matrix and its determinant associated with the map 𝐹_𝐾. Moreover, to each 𝐾∈𝒯, we associate the outward unit normal 𝑛_𝐾.

Usually, for DG methods, 𝐸 and 𝐻 fields are approximated by polynomials on each cell. In our case, we approximate by polynomials the fields 𝐷𝐹^∗_𝐾𝐸∘𝐹_𝐾 and 𝐷𝐹^∗_𝐾𝐻∘𝐹_𝐾 on 𝐾. It is not a strange choice since the Jacobian matrix is the essential ingredient to build a conform H-curl approximation [7, 8]. As we will see later, this will imply interesting properties for memory storage.

Finally, we consider the following semidiscrete DG method.

Find (𝐸_ℎ(⋅,𝑡),𝐻_ℎ(⋅,𝑡))∈𝐕_𝑟×𝐕_𝑟 such that for all 𝐾∈𝒯 and for all 𝜓,𝜙∈𝐕_𝑟𝐾𝜀𝜕_𝑡𝐸_ℎ⋅𝜓𝑑𝑥+_𝐾𝜎𝐸_ℎ=⋅𝜓𝑑𝑥_𝐾∇×𝐻_ℎ⋅𝜓𝑑𝑥+_𝜕𝐾𝛼^𝐾_𝜕𝐾𝑛_𝐾×(𝐸_ℎ×𝑛_𝐾)^𝐾_𝜕𝐾+𝛽^𝐾_𝜕𝐾𝐻_ℎ×𝑛_𝐾^𝐾_𝜕𝐾⋅𝜓𝑑𝑠,_𝐾𝜇𝜕_𝑡𝐻_ℎ⋅𝜙𝑑𝑥=−_𝐾∇×𝐸_ℎ+⋅𝜙𝑑𝑥_𝜕𝐾𝛾^𝐾_𝜕𝐾𝐸_ℎ×𝑛_𝐾^𝐾_𝜕𝐾+𝛿^𝐾_𝜕𝐾𝑛_𝐾×𝐻_ℎ×𝑛_𝐾^𝐾_𝜕𝐾⋅𝜙𝑑𝑠, (3) where 𝛼^𝐾_𝜕𝐾, 𝛽^𝐾_𝜕𝐾, 𝛾^𝐾_𝜕𝐾, 𝛿^𝐾_𝜕𝐾 are parameters constant per face and ⟦𝑣⟧^𝐾_Γ=(𝑣|_𝐾′)_Γ−(𝑣|_𝐾)_Γ is the jump across the boundary Γ=𝐾^′∩𝐾. When Γ is a boundary face (i.e., Γ=𝐾∩𝜕Ω) then 𝐾^′ does not exist and we define ⟦𝑣⟧^𝐾_Γ=−(𝑣|_𝐾)_Γ.

The coefficients 𝛼^𝐾_𝜕𝐾,𝛽^𝐾_𝜕𝐾,𝛾^𝐾_𝜕𝐾,𝛿^𝐾_𝜕𝐾 are chosen such that (1) and (3) are equivalent problems (in the continuous sense) and to ensure a conservative formulation. Then, we obtain(i)for all Γ=𝐾∩𝐾^′,𝛼^𝐾_Γ=𝛼^𝐾′_Γ=0,𝛿^𝐾_Γ=𝛿^𝐾′_Γ=0,𝛽^𝐾_Γ=𝛽^𝐾′_Γ=−1/2,𝛾^𝐾_Γ=𝛾^𝐾′_Γ=1/2,(ii)for all Γ=𝐾∩𝜕Ω,𝛽^𝐾_Γ=0=𝛿^𝐾_Γ=𝛼^𝐾_Γ,𝛾^𝐾_Γ=1.

For the time discretization, as for the FDTD method [9, 10], we use a Leap-Frog numerical scheme where the electric fields are evaluated at the time 𝑛Δ𝑡 and the magnetic fields at the time (𝑛+1/2)Δ𝑡, with Δ𝑡 defining the time step and 𝑛 the current iteration in time.

In order to define a set ℬ of basis functions of 𝐕_𝑟, we first define a set of basis functions of 𝐾. Let 𝐱_𝑖𝑗𝑘=(𝑥_𝑖,𝑦_𝑗,𝑧_𝑘), 1⩽𝑖,𝑗,𝑘⩽𝑟+1 be a set of points of 𝐾, where 𝑥_𝑖,𝑦_𝑗, and 𝑧_𝑘 are Gauss quadrature points on [0,1]. At the point 𝐱_𝑖𝑗𝑘, we define on 𝐾 three basis functions 𝜑^𝑙_𝑖𝑗𝑘(𝑥,𝑦,𝑧)=𝐿_𝑖(𝑥)𝐿_𝑗(𝑦)𝐿_𝑘(𝑧)𝑒^𝑙, where 𝐿_𝑖(𝑥)=∏^𝑟+1_{𝑚=1,𝑚≠𝑖}(𝑥−𝑥_𝑚)/(∏^𝑟+1_{𝑚=1,𝑚≠𝑖}(𝑥_𝑖−𝑥_𝑚)) is the Lagrange interpolation polynomial and (𝑒^𝑙)_𝑙=1,2,3 denotes the classical cartesian base. On 𝐾∈𝒯, the corresponding basis functions are defined by 𝜑^𝑙,𝐾_𝑖𝑗𝑘∘𝐹_𝐾(𝐱)=(𝐷𝐹^∗_𝐾)⁻¹𝜑^𝑙_𝑖𝑗𝑘(𝐱), where 𝐱=(𝑥,𝑦,𝑧). Finally, the set of basis function of 𝐕_𝑟 is ℬ=𝜑^𝑙,𝐾_𝑖𝑗𝑘∶𝐾∈𝒯,𝑙=1,2,3,𝑖,𝑗,𝑘=1,…,𝑟+1 (4) and the dimension of 𝐕_𝑟 is 3(𝑟+1)³𝑁 (where 𝑁 is the number of cells). So, each 𝑈(⋅,𝑡)∈𝐕_𝑟 can be written in this way: for all 𝐾∈𝒯𝑈_|𝐾(𝐱,𝑡)=³_𝑙=1^𝑟+1_{𝑖𝑗𝑘=1}𝑈(𝑡)^𝑙,𝐾_𝑖𝑗𝑘𝜑^𝑙,𝐾_𝑖𝑗𝑘(𝐱), (5) where 𝑈(𝑡)^𝑙,𝐾_𝑖𝑗𝑘∈ℝ is the degree of freedom associated to the basis function 𝜑^𝑙,𝐾_𝑖𝑗𝑘.

In the sequel, we say the DG method is a 𝑄_𝑟 approximation when we choose a spatial approximation order of 𝑟.

Thanks to the chosen approximation space, we have for all 𝐾∈𝒯,  𝑙,𝑠=1,2,3 and 𝑖,𝑗,𝑘,𝑚,𝑛,𝑝=1,…,𝑟+1, _𝐾∇×𝜑^l,𝐾_𝑖𝑗𝑘⋅𝜑^𝑠,𝐾_𝑚𝑛𝑝𝑑𝑥=sign𝐽_𝐾𝐾∇×𝜑^𝑙_𝑖𝑗𝑘⋅𝜑^𝑠_𝑚𝑛𝑝𝑑𝑥,_𝜕𝐾𝜑^𝑙,𝐾_𝑖𝑗𝑘×𝑛_𝐾⋅𝜑^𝑠,𝐾_𝑚𝑛𝑝𝑑𝑠=sign𝐽_𝐾𝜕𝐾𝜑^𝑙_𝑖𝑗𝑘×𝑛⋅𝜑^𝑠_𝑚𝑛𝑝𝑑𝑠,̂ (6) where 𝑛 is the outward unit normal to 𝐾.

By using (6) and a Gauss quadrature rule to evaluate integrals, we obtain the following:(i)for mass matrices, ((𝑀)_𝑝𝑙 denotes the (𝑝,𝑙) component of the matrix 𝑀): _𝐾𝜀𝑈⋅𝜑^𝑙,𝐾_𝑖𝑗𝑘𝑑𝑥=𝜔³_{𝑖𝑗𝑘𝑝=1}𝑈^𝑝,𝐾_𝑖𝑗𝑘𝐽_𝐾𝐷𝐹_𝐾⁻¹𝜀∘𝐹_𝐾𝐷𝐹_𝐾^∗−1_𝑝𝑙𝑥_𝑖𝑗𝑘, (7)(ii)for stiffness matrices, ((𝑀)_𝑙 denotes the component 𝑙 of vector 𝑀): _𝐾∇×𝑈⋅𝜑^𝑙,𝐾_𝑖𝑗𝑘=sign𝐽_𝐾^𝑟+13_{𝑚𝑛𝑞=1𝑝=1}𝜔_𝑖𝑗𝑘𝑈^𝑙,𝐾_𝑖𝑗𝑘∇×𝜑^𝑚_𝑖𝑗𝑘𝑥_{𝑖𝑗𝑘𝑙}, (8)(iii)for jump matrices _𝜕𝐾𝐻×𝑛_𝐾𝜕𝐾⋅𝜑^𝑙,𝐾_𝑖𝑗𝑘𝐽𝑑𝑠=−sign_𝐾′𝜕𝐾𝐷𝐹_𝐾⁻¹∘𝑆_𝐾′𝐾𝐷𝐹_𝐾′×𝐻_𝐾′×𝑛⋅𝜑^𝑙_𝑖𝑗𝑘∘𝑆_𝐾′𝐾𝑑𝑠̂+sign𝐽_𝐾𝜕𝐾𝐻_𝐾×𝑛⋅𝜑^𝑙_𝑖𝑗𝑘𝑑𝑠.̂ (9)In the above expressions 𝜔_𝑖𝑗𝑘 is the quadrature weight at point 𝑥_𝑖𝑗𝑘 and 𝑆_𝐾′𝐾=𝐹_𝐾⁻¹∘𝐹_𝐾′. 𝐷𝐹_𝐾⁻¹∘𝑆_𝐾′𝐾𝐷𝐹_𝐾′ is a permutation matrix constant per face [1].

The DG formulation (3) finally leads to the semidiscrete formulation: 𝑀_𝜀𝜕_𝑡𝐸+𝑀_𝜎𝐸=𝑅𝐻−𝑆^𝑖𝑀𝐻,_𝜇𝜕_𝑡𝐻=−𝑅𝐸+𝑆^𝑖𝐸+𝑆^𝑏𝐸, (10) where 𝑀_𝜀,𝑀_𝜇, and 𝑀_𝜎 are 3×3 block-diagonal matrices, 𝑅, the stiffness matrix, and 𝑆^𝑖 and 𝑆^𝑏 are jump matrices. Thanks to the choice of the space approximation and of the basis functions, the mass matrix is given, for all spatial orders, by a 3×3 block-diagonal matrix which needs to be stored for each cell 𝐾. Stiffness and jump matrices just require to store the sign of the Jacobians 𝐽_𝐾 and some computations made on the reference element 𝐾.

Then, the important advantage of this DG method is to give, regardless of the space approximation order, a very low memory storage and a small cost of computation to evaluate the matrices of the numerical scheme. This allows us to use meshes with a small number of cells and a high order spatial approximation to obtain very accurate solutions. Consequently, to obtain an accurate solution, the memory and CPU costs needed to solve a problem with this approach are lower than those for some other methods, as, for example, FDTD, based on a spatial approximation of order 2 and using more refined meshes. In fact, high order spatial approximation reduces the dispersive and dissipative errors of the scheme and improve the accuracy of the solution. The use of non-cubic cells, like for the FDTD method, in the meshes in the method permits also to evaluate correctly the fields near the structures. Therefore the proposed method is well adapted not only to Electromagnetic coupling problems for which it is essential to know this kind of values, but also to cavity problems where the dispersive and dissipative errors may not be neglected.

2.2. Local Time-Stepping Method

The previous DG numerical scheme studied is based for the time domain on an explicit Leap-Frog discretization [9]. Nevertheless, when dealing with unstructured meshes, we can obtain significant cell size disparities and, to ensure stability, the global time step is constrained by the smallest cell of the mesh. This implies an important increase of the computational time. To avoid this problem, a local time-stepping strategy has been proposed for the Leap-Frog time discretization [2]. This method is based on a decomposition of the cells of the mesh in a set of 𝑁 classes of cells where each class 𝑖 has a time step given by 3^𝑁−𝑖Δ𝑡_min𝑟, with Δ𝑡_min defining the smallest time step on all the classes.

For example, when the problem has 3 classes, the different operations that proceed in a step of the local time-stepping Leap-Frog method are given in Figure 1. We consider that the time step for the method is given by the largest value on all the classes.

Figure 1 

One Step of the multiclass Leap-Frog method.

This local time-stepping strategy can be easily implemented as a recursive process.

3. Mixed-Spatial Order Approximation Strategy

In the mixed-spatial order strategy proposed, we search to evaluate a solution of a problem by considering a given accuracy. Generally, to treat a problem, the mesh of the geometry is defined by taking into account the largest frequency of investigation contained in the spectrum of the source. Indeed, the size of the cells is smaller than 𝜆_min/𝑛 where 𝜆_min defines the smallest wavelength of the source and 𝑛 is an integer which depends on the wanted accuracy. For example, in an FDTD approach, generally, people use 𝑛=10. In other terms, the length 𝜆_min/𝑛 corresponds to the distance between 2 degrees of freedom in the numerical method. In the DG method, the number of degrees of freedom in each cell depends on the order of the considered spatial approximation. In particular, for a given mesh and a given spatial order 𝑟, to keep a discretization of 𝜆_min/𝑛 between each degree of freedom, the ratio of the size of the cell by the smallest wavelength will be given by 𝑟/𝑛. For example, by considering 𝑛=10, for a 𝑄1 approximation we obtain 0.1 wavelength for the ratio of the size of the cells, for 𝑄2 we obtain 0.2, and for 𝑄3 we obtain 0.3. Then, by knowing the size of a cell and the smallest wavelength of the source, we can determine the spatial degree necessary on each cell to maintain a global discretization of 𝜆_min/𝑛 between each degree of freedom. The mixed-spatial order strategy proposed is the following: (i)define the smallest wavelength 𝜆_min by considering the spectrum of the source; (ii)define an equivalent geometric length ℎ on each cell; (iii)allocate a spatial order 𝑟=(ℎ/𝜆_min)𝑛 on each cell by considering ℎ, 𝜆_min, and a given 𝑛.

The difficulty in this strategy is to determine a value of ℎ for each cell equivalent to its size. Several possibilities have been tried and the most interesting is a value which corresponds to the largest diagonal length of the element. By considering this value, we do not overestimate or underestimate the spatial discretization in any direction of the element. The example given in Figure 2 shows the gain that we can expect by applying a mixed-spatial order strategy. In particular, we can see that the time step of the method using mixed-spatial order is higher than that without using it. This implies a smaller number of time iterations to describe a given range in time and then a gain in computational time. We note also a decrease in the number of degrees of freedom in the problem which leads to a lower required memory. However to confirm the interest of this approach, it is also important to compare the accuracy of the solutions obtained. This point has been developed in Section 3.3.

Figure 2 

Potential gain in memory storage and CPU time by using a mixed-spatial order strategy in the DG method.

3.1. Jump Matrix

To implement the mixed-spatial order strategy proposed, we only need to modify the calculation of the “crossed-jump” terms in the DG scheme given in Section 2. Indeed, as it has been described in the second section, the evaluation of the mass and the stiffness matrices is local at each cell level. Then, the mixed-order strategy proposed does not modify the computation of these matrices. The jump terms are split into 2 terms where the calculation of one of them is unchanged because the values considered to compute it depend only on one element. In the calculation of the second term, the values taken into account depend on two elements and so the evaluation of this term must be changed. As seen in the first section, we have to use a Gauss quadrature formula to evaluate the integral terms. To ensure accurate evaluation of the “crossed jump” terms on each adjacent cell to the face considered, we use, in the quadrature formula, a number of point’s corresponding to the largest spatial order of the two cells. Let 𝜑^𝑙0_𝑖^,𝐾_0𝑗0𝑘0 be a given basis function on 𝐾; by considering two contiguous cells 𝐾 and 𝐾^′ with different spatial approximated orders 𝑟 and 𝑟^′, 𝑟>𝑟^′, the jump term for the equation of the magnetic field in the system (3) on the cell 𝐾 is split into a centered term 𝑆_𝑐𝐾=∫_Γ⟦𝐸×𝑛_𝐾⟧⋅𝜑^𝑙0_𝑖^,𝐾_0𝑗0𝑘0𝑑𝑠 and a dissipative term 𝑆_𝑑𝐾=∫_Γ⟦𝑛_𝐾×(𝐻×𝑛_𝐾)⟧⋅𝜑^𝑙0_𝑖^,𝐾_0𝑗0𝑘0𝑑𝑠, where Γ=𝐾∩𝐾^′. For the evaluation of 𝑆_𝑐𝐾=∫_Γ(𝐸_𝐾′×𝑛_𝐾)⋅𝜑^𝑙0_𝑖^,𝐾_{0,𝑗0,𝑘0}𝑑𝑠−∫_Γ(𝐸_𝐾×𝑛_𝐾)⋅𝜑^𝑙0_𝑖^,𝐾_0𝑗0𝑘0𝑑𝑠, if we consider 𝑛_𝐾 as the unit outward normal to the cell 𝐾, we have ∀𝑥∈Γ,𝑛_𝐾∘𝐹_𝐾(𝑥)=𝐷𝐹_𝐾^∗−1(𝑥)𝑛𝐷𝐹_𝐾^∗−1(𝑥)𝑛,𝑑𝑠=𝐷𝐹_𝐾^∗−1(𝑥)𝑛𝐽_𝐾(𝑥), (11) where 𝑥=𝐹_𝐾(𝑥). The term 𝐸_𝐾(𝑥)×𝑛_𝐾 can be written as 𝐷𝐹_𝐾^∗−1(𝑥)𝐸(𝑥)×𝐷𝐹_𝐾^∗−1(𝑥)𝑛𝐷𝐹_𝐾^∗−1(𝑥)𝑛, (12) where 𝐸(𝑥)=³_𝑙=1^𝑟+13_{𝑖𝑗𝑘=1}𝐸^𝑙,𝐾_𝑖𝑗𝑘𝜑^𝑙_𝑖𝑗𝑘(𝑥). (13) This term is also equal to 𝐷𝐹_𝐾(𝑥)𝐽_𝐾(𝑥)𝐸(𝑥)×𝑛. (14) Then we obtain _Γ𝐸_𝐾×𝑛_𝐾⋅𝜑^𝑙0_𝑖^,𝐾_0𝑗0𝑘0=Γ𝐽𝑑𝑠sign_𝐾(𝑥)𝐸(𝑥)×𝑛⋅𝜑^𝑙0_{𝑖0𝑗0𝑘0}(𝑥)𝑑𝑠.̂ (15) Concerning the term ∫_Γ(𝐸_𝐾′×𝑛_𝐾)⋅𝜑^𝑙0_𝑖^,𝐾_0𝑗0𝑘0𝑑𝑠, we take 𝑆_𝐾𝐾′=𝐹_𝐾⁻¹_′∘𝐹_𝐾 with 𝐹_𝐾′ being the trilinear transformation between the cell 𝐾^′ and the reference element 𝐾. Then, by using 𝑛_𝐾=−𝑛_𝐾′, the previous integral term can be written as −Γ𝐷𝐹_𝐾′∘𝑆_𝐾𝐾′(𝑥)𝐽_𝐾′∘𝑆_𝐾𝐾′(𝑥)𝐸_𝐾′×𝑛^′𝐷𝐹_𝐾^∗−1_′(𝑥)𝑛^′⋅𝐷𝐹_𝐾^∗−1𝜑^𝑙0_{𝑖0𝑗0𝑘0}(𝑥)𝐷𝐹_𝐾^∗−1(𝑥)𝑛𝐽_𝐾(𝑥)𝑑𝑠̂ (16) with 𝑛^′ defined as 𝑛_𝐾′=𝐷𝐹_𝐾^∗−1_′∘𝑆_𝐾𝐾′(𝑥)𝑛^′𝐷𝐹_𝐾^∗−1_′∘𝑆_𝐾𝐾′(𝑥)𝑛^′. (17) By using the relation 𝐷𝐹_𝐾^∗−1(𝑥)𝑛𝐽_𝐾(𝑥=)𝐷𝐹_𝐾^∗−1_′∘𝑆_𝐾𝐾′(𝑥)𝑛^′𝐽_𝐾′∘𝑆_𝐾𝐾′(𝑥) (18) we obtain for the previous term −Γsign𝐽_𝐾∘𝑆_𝐾𝐾′(𝑥)𝐷𝐹_𝐾⁻¹(𝑥)𝐷𝐹_𝐾′∘𝑆_𝐾𝐾′(𝑥)𝐸_𝐾′(𝑥)×𝑛^′⋅𝜑^𝑙0_{𝑖0,𝑗0,𝑘0}(𝑥)𝑑𝑠.̂ (19) Finally, for the term 𝑆_𝑐𝐾, we have𝑆_𝑐𝐾Γ=−sign𝐽_𝐾∘𝑆_𝐾𝐾′(𝑥)𝐷𝐹_𝐾⁻¹(𝑥)𝐷𝐹_𝐾′∘𝑆_𝐾𝐾′(𝑥)𝐸_𝐾′(𝑥)×𝑛^′⋅𝜑^𝑙0_{𝑖0,𝑗0,𝑘0}(𝑥)𝑑𝑠̂−Γsign𝐽_𝐾(𝑥)𝐸(𝑥)×𝑛⋅𝜑^𝑙0_{𝑖0𝑗0𝑘0}(𝑥)𝑑𝑠.̂ (20)

For the dissipative term we have 𝑆_𝑑𝐾=_Γ𝑛_𝐾×𝐻_𝐾′×𝑛_𝐾⋅𝜑^𝑙0_𝑖^,𝐾_0𝑗0𝑘0−𝑑𝑠_Γ𝑛_𝐾×𝐻_𝐾×𝑛_𝐾⋅𝜑^𝑙0_𝑖^,𝐾_0𝑗0𝑘0𝑑𝑠. (21) By using 𝐻(𝑥)=∑³_𝑙=1∑^(𝑟+1)3_{𝑖𝑗𝑘=1}𝐻^𝑙,𝐾_𝑖𝑗𝑘𝜑^𝑙_𝑖𝑗𝑘(𝑥), the second term of (21) can be written as _Γ𝐻_𝐾×𝑛_𝐾⋅𝜑^𝑙0_𝑖^,𝐾_0𝑗0𝑘0×𝑛_𝐾Γ𝑑𝑠=𝐷𝐹_𝐾(𝑥)𝐽_𝐾(𝑥)𝐻_𝐾(𝑥)×𝑛𝐷𝐹_𝐾^∗−1(𝑥)𝑛⋅𝐷𝐹_𝐾(𝑥)𝐽_𝐾(𝑥)𝜑^𝑙0_{𝑖0𝑗0𝑘0}×𝑛𝐷𝐹_𝐾^∗−1(𝑥)𝑛𝐽_𝐾(𝑥𝑑𝑠̂=Γ𝐽)sign_𝐾(𝑥)𝐷𝐹^∗_𝐾(𝑥)𝐷𝐹_𝐾(𝑥)𝐽_𝐾(𝑥)𝐷𝐹_𝐾^∗−1(𝑥)𝑛×𝐻_𝐾×𝑛⋅𝜑^𝑙0_{𝑖0𝑗0𝑘0}×𝑛𝑑𝑠.̂ (22) In the same way, the first term of (21) becomes Γ𝐷𝐹_𝐾(𝑥)𝐷𝐹_𝐾′∘𝑆_𝐾𝐾′(𝑥)𝐽_𝐾′∘𝑆_𝐾𝐾′(𝑥)𝐷𝐹_𝐾^∗−1_′∘𝑆_𝐾𝐾′(𝑥)𝑛^′×𝐻_𝐾′×𝑛^′⋅𝜑^𝑙0_{𝑖0𝑗0𝑘0}×𝑛𝑑𝑠.̂ (23) Finally, we obtain 𝑆_𝑑𝐾=Γ𝐷𝐹_𝐾(𝑥)𝐷𝐹_𝐾′∘𝑆_𝐾𝐾′(𝑥)𝐽_𝐾′∘𝑆_𝐾𝐾′(𝑥)𝐷𝐹_𝐾^∗−1_′∘𝑆_𝐾𝐾′(𝑥)𝑛^′×𝐻_𝐾′×𝑛^′⋅𝜑^𝑙0_{𝑖0𝑗0𝑘0}×𝑛𝑑𝑠̂Γ−sign𝐽_𝐾(𝑥)𝐷𝐹^∗_𝐾(𝑥)𝐷𝐹_𝐾(𝑥)𝐽_𝐾(𝑥)𝐷𝐹_𝐾^∗−1(𝑥)𝑛×𝐻_𝐾×𝑛⋅𝜑^𝑙0_{𝑖0𝑗0𝑘0}×𝑛𝑑𝑠.̂ (24)

Similarly, the jump term on the cell 𝐾^′ is given by a centered term 𝑆_𝑐𝐾′=∫_Γ⟦𝐸×𝑛^′⟧⋅𝜑^′𝑑𝑠 and a dissipative term 𝑆_𝑑𝐾′=∫_Γ⟦𝑛^′×(𝐻×𝑛^′)⟧⋅𝜑^′𝑑𝑠. By considering a basis function 𝜑^{𝑙0,𝐾′}_{𝑖0,𝑗0,𝑘0} on 𝐾^′ and a similar development as for the jump terms of the cell 𝐾, we obtain 𝑆_𝑐𝐾′Γ=−sign𝐽_𝐾′∘𝑆_𝐾′𝐾(𝑥)𝐷𝐹_𝐾⁻¹_′(𝑥)𝐷𝐹_𝐾∘𝑆_𝐾′𝐾(𝑥)𝐸_𝐾(𝑥)×𝑛^′⋅𝜑^𝑙0_{𝑖0,𝑗0,𝑘0}(𝑥)𝑑𝑠̂−Γsign𝐽_𝐾′(𝑥)𝐸_𝐾′(𝑥)×𝑛⋅𝜑^𝑙0_{𝑖0𝑗0𝑘0}(𝑥)𝑑̂𝑆𝑠,_𝑑𝐾′=Γ𝐷𝐹_𝐾′(𝑥)𝐷𝐹_𝐾∘𝑆_𝐾′𝐾(𝑥)𝐽_𝐾∘𝑆_𝐾′𝐾(𝑥)𝐷𝐹_𝐾^∗−1∘𝑆_𝐾′𝐾(𝑥)𝑛^′×𝐻_𝐾×𝑛^′⋅𝜑^𝑙0_{𝑖0𝑗0𝑘0}×𝑛𝑑𝑠̂−Γsign𝐽_𝐾′(𝑥)𝐷𝐹^∗_𝐾′(𝑥)𝐷𝐹_𝐾′(𝑥)𝐽_𝐾′(𝑥)𝐷𝐹_𝐾^∗−1_′(𝑥)𝑛×𝐻_𝐾′×𝑛⋅𝜑^𝑙0_{𝑖0𝑗0𝑘0}×𝑛𝑑𝑠,̂ (25) where Γ=𝐹_𝐾′Γ=𝐹_𝐾Γ^′,𝑆_𝐾′𝐾=𝐹_𝐾^∗−1∘𝐹_𝐾′,𝑥=𝐹_𝐾⁻¹_′(𝑥),𝑛_𝐾′∘𝐹_𝐾′(𝑥)=𝐷𝐹_𝐾^∗−1_′(𝑥)𝑛𝐷𝐹_𝐾^∗−1_′(𝑥)𝑛,𝑛_𝐾∘𝐹_𝐾′(𝑥)=𝐷𝐹_𝐾^∗−1∘𝑆_𝐾′𝐾(𝑥)𝑛^′𝐷𝐹_𝐾^∗−1∘𝑆_𝐾′𝐾(𝑥)𝑛^′. (26)

By using a Gauss quadrature formula of order equal to the largest spatial approximation order of 𝐾 and 𝐾^′ to evaluate the integrals, our choice of basis functions leads to having only few nonnull interactions. Figure 3 in 2D case, gives the nonnull contributions of the different basis functions considering a same spatial order approximation for both adjacent elements to the face. In this case, for each point of quadrature, only a line of basis functions is used (this is also true in 3D). This is different when the adjacent cells to the considered face have not the same spatial order approximation. In this case, to compute the “crossed-jump” terms, we need to use all the basis functions for each point of quadrature (2D case as 3D case). Figures 4 and 5 show examples of the different interactions taken into account to compute these terms. We show in these figures that the calculation processes are not symmetric.

Figure 3 

Jump terms by considering 𝑟=𝑟′=3.

Figure 4 

Jump terms associated to 𝐾, by considering 𝑟=3 and 𝑟′=2.

Figure 5 

Jump terms associated to 𝐾′, by considering 𝑟=3 and 𝑟′=2.

3.2. Stability

In this part, we are reinvestigating the stability results already obtained in [1, 2] for a nondissipative and a dissipative DG formulation with a same spatial order of approximation for all the cells in the computational domain. These studies are based on the analysis of a discrete energy. For example, in the case of the nondissipative scheme, one defines the discrete quantityℰ^𝑛_ℎ=_{𝐾∈𝒯ℎ}^𝐺_𝐾𝜀_𝐾𝐸^𝑛_ℎ,𝐾⋅𝐸^𝑛_ℎ,𝐾𝑑𝑥+^𝐺_𝐾𝜇_𝐾𝐻^𝑛+1/2_ℎ,𝐾⋅𝐻^𝑛−1/2_ℎ,𝐾𝑑𝑥 (27) (where ∫^𝐺_𝐾 specifies that one that uses a Gauss quadrature rule to compute the integral) which is conserved during the time, that is, ℰ_ℎ^𝑛+1=ℰ^𝑛_ℎ, and which becomes a discrete energy under an explicit CFL condition.

In the case of the mixed-spatial order strategy proposed in this paper, this result still holds.

Proposition 1. For the multiorder nondissipative scheme, the discrete quantity ℰ^𝑛_ℎ is conserved during the discrete time: ∀𝑛≥0,ℰ_ℎ^𝑛+1−ℰ^𝑛_ℎ=0. (28)

Proof. This proof is classical and does not raise any difficulty. This is why we give only the key points of the proof. First, by taking the test functions 𝜓=𝐸_ℎ^𝑛+1+𝐸^𝑛_ℎ and 𝜙=𝐻_ℎ^𝑛+1/2 in the nondissipative fully discrete scheme, one can easily obtain [5] ℰ_ℎ^𝑛+1−ℰ^𝑛_ℎ=_{𝐾∈𝒯𝐾}∇×𝐻_𝐾^𝑛+1/2⋅𝐸_𝐾^𝑛+1+𝐸^𝑛_𝐾−𝑑𝑥_𝐾∇×𝐸_𝐾^𝑛+1+𝐸^𝑛_𝐾⋅𝐻_𝐾^𝑛+1/2+𝑑𝑥_𝐺𝜕𝐾𝛽𝐻_ℎ^𝑛+1/2×𝑛_𝐾^𝐾_𝜕𝐾⋅𝐸_𝐾^𝑛+1+𝐸^𝑛_𝐾+^𝐺_𝜕𝐾𝛾𝐸_ℎ^𝑛+1+𝐸^𝑛_ℎ×𝑛_𝐾^𝐾_𝜕𝐾⋅𝐻_𝐾^𝑛+1/2=_{𝐾∈𝒯𝜕𝐾}𝑛_𝐾×𝐻_𝐾^𝑛+1/2⋅𝐸_𝐾^𝑛+1+𝐸^𝑛_𝐾+𝑑𝜎_𝐺𝜕𝐾𝛽𝐻_ℎ^𝑛+1/2×𝑛_𝐾^𝐾_𝜕𝐾⋅𝐸_𝐾^𝑛+1+𝐸^𝑛_𝐾+^𝐺_𝜕𝐾𝛾𝐸_ℎ^𝑛+1+𝐸^𝑛_ℎ×𝑛_𝐾^𝐾_𝜕𝐾⋅𝐻_𝐾^𝑛+1/2. (29) Now, we examine a face Γ=𝐾₊∩𝐾₋. We obtain the terms 𝑇=_Γ𝑛₊×𝐻₊^𝑛+1/2⋅𝐸₊^𝑛+1+𝐸^𝑛₊+𝑛₋×𝐻₋^𝑛+1/2⋅𝐸₋^𝑛+1+𝐸^𝑛₋+𝑑𝜎_𝐺+Γ𝛽𝐻₋^𝑛+1/2×𝑛₊−𝐻₊^𝑛+1/2×𝑛₊⋅𝐸₊^𝑛+1+𝐸^𝑛₊+^𝐺−_Γ𝛽𝐻₊^𝑛+1/2×𝑛₋−𝐻₋^𝑛+1/2×𝑛₋⋅𝐸₋^𝑛+1+𝐸^𝑛₋+^𝐺+_Γ𝛾𝐸₋^𝑛+1+𝐸^𝑛₋×𝑛₊−𝐸₊^𝑛+1+𝐸^𝑛₊×𝑛₊⋅𝐻₊^𝑛+1/2+^𝐺−_Γ𝛾𝐸₊^𝑛+1+𝐸^𝑛₊×𝑛₋−𝐸₋^𝑛+1+𝐸^𝑛₋×𝑛₋⋅𝐻₋^𝑛+1/2. (30) We can decompose 𝑇 under the form 𝑇1^±+𝑇2+𝑇3, where 𝑇1^±=_Γ𝑛_±×𝐻_±^𝑛+1/2⋅𝐸_±^𝑛+1+𝐸^𝑛_±−𝑑𝜎_𝐺±Γ𝛽𝐻_±^𝑛+1/2×𝑛_±⋅𝐸_±^𝑛+1+𝐸^𝑛_±−^𝐺±_Γ𝛾𝐸_±^𝑛+1+𝐸^𝑛_±×𝑛_±⋅𝐻_±^𝑛+1/2,𝑇2=^𝐺+_Γ𝛽𝐻₋^𝑛+1/2×𝑛₊⋅𝐸₊^𝑛+1+𝐸^𝑛₊−^𝐺−_Γ𝛾𝐸₊^𝑛+1+𝐸^𝑛₊×𝑛₋⋅𝐻₋^𝑛+1/2,𝑇3=^𝐺−_Γ𝛽𝐻₊^𝑛+1/2×𝑛₋⋅𝐸₋^𝑛+1+𝐸^𝑛₋−^𝐺+_Γ𝛾𝐸₋^𝑛+1+𝐸^𝑛₋×𝑛₊⋅𝐻₊^𝑛+1/2. (31)
So, if one takes for 𝑇1^± the Gauss quadrature rule of order 𝑟_± and for 𝑇2,𝑇3 the same quadrature rule for 𝐺₊ and 𝐺₋ (not necessarily the one corresponding to the highest spatial order), one obtains 𝑇=0.

So, one can give the stability result [1].

Theorem 2. Let ℛ_𝐾,𝒟_𝐾,ℬ_𝐾 be the 3(𝑟_𝐾+1)³×3(𝑟_𝐾+1)³ matrices defined by the following: for all 𝑙,𝑙^′∈{1,2,3} and for all 𝐼=(𝑖,𝑗,𝑘), 𝐼^′=(𝑖^′,𝑗^′,𝑘^′)∈{1,…,𝑟_𝐾+1}³ℛ_𝐾(𝑙,𝐼),𝑙^′,𝐼^′=_𝜕𝐾∇×𝜑^𝑙_𝑖𝑗𝑘⋅∇×𝜑^𝑙′_{𝑖′𝑗′𝑘′},𝒟_𝐾(𝑙,𝐼),𝑙^′,𝐼^′=𝛿_𝑙𝑙′𝛿_𝐼𝐼′𝜔_𝑖𝑗𝑘,ℬ_𝐾(𝑙,𝐼),𝑙^′,𝐼^′=_𝜕𝐾𝜑^𝑙_𝑖𝑗𝑘×𝑛⋅𝜑^𝑙′_{𝑖′𝑗′𝑘′}×𝑛𝑑𝑠.̂ (32) The condition given by Δ𝑡Λ_𝐾<2𝑐_𝐾1×𝜆_max𝒟_𝐾^−1/2ℛ_𝐾𝒟_𝐾^−1/2+12max_{1≤𝑖≤𝑛𝑏𝑓𝑖𝐾}𝜇_𝐾𝜇_{𝑉(𝑖,𝐾)},𝜀_𝐾𝜀_{𝑉(𝑖,𝐾)}×𝜆_max𝒟_𝐾^−1/2ℬ_𝐾𝒟_𝐾^−1/2−1 (33) for all 𝐾 is sufficient to ensure the stability of the DG scheme; where Λ_𝐾=min_{1≤𝑖,𝑗,𝑘≤𝑟𝐾+1}𝐽_𝐾𝐱_𝑖𝑗𝑘𝜆_max𝐷𝐹^∗_𝐾𝐷𝐹_𝐾𝐱_𝑖𝑗𝑘 (34) (in the case of a uniform cartesian grid whose spatial step is ℎ, Λ_𝐾=ℎ), where 𝑛𝑏𝑓𝑖_𝐾 is the number of internal faces of 𝐾 and 𝑉(𝑖,𝐾) is the the neighbor of 𝐾.

Proof. The proof of the theorem is the same as the one given in the paper [1]. We must only take into account a different spatial order for each cell. However this does not imply any modification in the proof.

Remark 3. In the case of the dissipative scheme, the stability can be also proved by using a modified discrete energy analysis as defined in [2].

3.3. Numerical Results

In this section, we present three examples to quantify the advantages of the proposed mixed-spatial solution strategy. The first example consists in evaluating the field scattered by a 1𝑚×1𝑚 perfectly metallic surface illuminated with a plane wave. This evaluation has been done for different meshes where the ratio between the smallest and the largest cell is progressively increased in order to generate several spatial approximation orders. Figures 6, 7, 8, and 9 show partial views of the meshes, where for the same accuracy of the solution we compare the CPU-time and the memory storage. The initial mesh (Mesh1) consists in an uniform cartesian mesh of the computational domain. The other meshes are obtained from the initial mesh by recursively cutting one cell in two (see Figures 6, 7, 8, and 9). In this example, by considering the spectrum of the source and the largest size of the cell, our discretization criterion (see Section 3) leads to a maximal spatial order 𝑟=7. So, we obtain a 𝑄₇ approximation for our meshes. Figures 10 and 11 show the CPU-time and the memory storage necessary for the DG method used with and without the mixed-spatial order strategy. We observe an important gain in CPU-time and in memory storage by using the mixed-spatial order approach. In particular we can note that our strategy implies a CPU time and a memory storage which are almost mesh independent.

Figure 6 

Partial view of Mesh1: a 3×3×3 uniform cartesian mesh.

Figure 7 

Partial view of Mesh2: mesh obtained by cutting a cell of Mesh1 in two.

Figure 8 

Partial view of Mesh3: mesh obtained by cutting a cell of Mesh2 in two.

Figure 9 

Partial view of Mesh4: mesh obtained by cutting a cell of Mesh3 in two.

Figure 10 

Evolution of the CPU-time for the different meshes used in the computation.

Figure 11 

Evolution of the memory storage for the different meshes used in the computation.

The second example consists in evaluating the propagation of a mode inside a cavity for a long time (5 μs). The considered cavity is a perfectly metallic cube where the length of the sides is taken to one meter. The studied mode is given by 𝐸_𝑥𝐸=0,_𝑦𝐸=0,_𝑧()(),𝐻=sin𝑚𝜋𝑥sin𝑛𝜋ycos𝜔𝑡_𝑥1=−𝜔𝜇₀()(),𝐻𝜋𝑛sin𝑚𝜋𝑥cos𝑛𝜋𝑦sin𝜔𝑡_𝑦=1𝜔𝜇₀()(),𝐻𝜋𝑚cos𝑚𝜋𝑥sin𝑛𝜋𝑦sin𝜔𝑡_𝑧=0, (35) where 𝜔=𝑐₀(𝑚𝜋)²+(𝑛𝜋)², 𝑚=1, and 𝑛=1.

In this example, we evaluate the field at the center of the cavity for a given nonstructured mesh and we compare the solutions given by the DG method using or not the mixed-spatial order strategy. For our mixed-spatial order strategy, we obtain two orders (3 and 4) which are assigned in the mesh as shown in Figure 12.

Figure 12 

Distributions for 𝑄3 and 𝑄4 cells projected on a (𝑥,𝑦) plane.

Figure 13 shows the analytical solution and the solutions obtained by the 𝑄3, 𝑄4, and the mixed 𝑄3-𝑄4 DG approximations. We can see that the 𝑄3-𝑄4 DG approximation leads to the same accurate solution of the 𝑄4 DG method but its numerical costs are far from smaller (see Table 1). Finally, we show in Figure 14 that the 𝑄3 solution is not sufficiently accurate whereas the mixed 𝑄3-𝑄4 approximation gives the best compromise between CPU-time, memory storage, and accuracy.

Table 1 

Time CPU and memory storage comparison between 𝑄3, 𝑄4, and mixed-spatial 𝑄3-𝑄4 approaches.

(a)

(b)

Figure 13 

Comparison between 𝑄3, 𝑄4, mixed-spatial order, and analytic solutions for the electric field 𝐸𝑧 taken at the center of the cavity. The ratio defines the percent of cells for a given order in the mesh. For example in mixed 𝑄3-𝑄4, we have 87.5% of 𝑄3 cells and 12.5% of 𝑄4 cells.

(a)

(b)

Figure 14 

Point to point error for the different computed solutions and L2 error.

In the last example we propose to evaluate the field scattered by an aircraft illuminated with a plane wave, by using the local time-stepping strategy given in Section 2.2, jointly with the mixed-spatial order strategy.

By considering the mesh (see Figure 15) and the spectrum of the source, the maximal order detected by our static adaptive strategy is equal to 2 in order to obtain 10 points per wavelength for the spatial discretization. So, theoretically, one must use a 𝑄2 approximation to ensure an accurate computed solution if we consider only one spatial order for the DG method. However, due to the difference of the cells sizes in the mesh, the 𝑄1 spatial approximation could be sufficient in some parts of the mesh. Then by using a mixed 𝑄1-𝑄2 spatial order approximation, we obtain for this example an important gain in CPU-time (factor 7) and in memory storage (see Table 2). Figures 16 and 17 show the electric and magnetic fields at the test-point 𝐴 (see Figure 15) obtained by 𝑄1, 𝑄2, and 𝑄1-𝑄2 DG approximation for the same mesh. We can see that the solutions obtained by the mixed-spatial order is similar to the 𝑄2 DG approximation, with a best cost in CPU-time and memory storage (see Table 2).

Table 2 

Results obtained for the aircraft example.

Figure 15 

Test-point evaluated on the aircraft.

(a)

(b)

Figure 16 

Comparison on test point 𝐴 for the electric field.

(a)

(b)

Figure 17 

Comparison on test point 𝐴 for the magnetic field.

4. Conclusion

In this paper we presented a p-strategy with a local time-stepping strategy for a particular Discontinuous Galerkin method to solve Maxwell's equations in time domain. We proposed a simple strategy based on the frequency spectrum of the source and on a characteristic size of the cells to evaluate the spatial approximation order to be applied on each cell of the mesh. We showed the validation and the advantages of this method on simple examples. The stability of the method has been also studied and a stability condition has been established. Finally, this method has been applied to study the electromagnetic fields scattered by a 3D object typical of aeronautic applications and give good results with a very interesting gain in CPU-time and memory storage. However, in the proposed strategy, the order is the same in all the directions of the elements and this can be a drawback for elements stretched in a particular direction. In this case, it may be more interesting to affect an order for each direction inside each element. This is one objective of our future works.

Acknowledgment

This work has been supported by the French Délégation Générale de l'Armement (DGA).