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

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.


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 timestepping 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 timestepping 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.

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.

Discontinuous Galerkin Method.
Let Ω be a bounded open subset of R 3 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, [ → R 3 × R 3 such that 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 T of hexahedral elements (  ) =1⋅⋅⋅ being a partition of Ω.We introduce the following space of approximation: where K = [0, 1] 3 is the unit cube or the reference element, for all  ∈ T,   : K →  denotes the trilinear mapping which associates the vertices of each element,   ( K) is the space of polynomials of degree at most equal to  ∈ N * in each variable on K, and   and   are, respectively, the Jacobian matrix and its determinant associated with the map   .Moreover, to each  ∈ T, 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 K.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 ( ℎ (⋅, ),  ℎ (⋅, )) ∈ V  × V  such that for all  ∈ T and for all ,  ∈ V where    ,   3) are equivalent problems (in the continuous sense) and to ensure a conservative formulation.Then, we obtain 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 B of basis functions of V  , we first define a set of basis functions of K. Let x = (x  , ŷ , ẑ ), 1 ⩽ , ,  ⩽  + 1 be a set of points of K, where x , ŷ , and ẑ are Gauss quadrature points on [0, 1].At the point x , we define on K three basis functions φ  ( x, ŷ, ẑ) =   ( x)  ( ŷ)  (ẑ)  , where is the Lagrange interpolation polynomial and (  ) =1,2,3 denotes the classical cartesian base.On  ∈ T, the corresponding basis functions are defined by  ,  ∘   (x) = ( *  ) −1 φ  (x), where x = (x, ŷ, ẑ).Finally, the set of basis function of V  is and the dimension of V  is 3(+1) 3  (where  is the number of cells).So, each (⋅, ) ∈ V  can be written in this way: for all  ∈ T where () ,  ∈ R 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  ∈ T, ,  = 1, 2, 3 and , , , , ,  = 1, . . .,  + 1, where n is the outward unit normal to K. 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 ): (ii) for stiffness matrices, (()  denotes the component  of vector ): (iii) for jump matrices In the above expressions   is the quadrature weight at point x and     =  −1  ∘    . −1  ∘        is a permutation matrix constant per face [1].
The DG formulation (3) finally leads to the semidiscrete formulation: 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 K.
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.

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.This local time-stepping strategy can be easily implemented as a recursive process.

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.

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 where Γ =  ∩   .For the evaluation of  where  =   ( x).The term   () ×   can be written as where This term is also equal to Then we obtain Concerning the term ∫ Γ (   ×   ) ⋅   0 ,  0  0  0 , we take    =  −1   ∘   with    being the trilinear transformation between the cell   and the reference element K.Then, by using   = −   , the previous integral term can be written as with   defined as By using the relation we obtain for the previous term Finally, for the term   , we have For the dissipative term we have By using Ĥ( x) = ∑ 3 =1 ∑ (+1) 3 =1  ,  φ  ( x), the second term of (21) can be written as In the same way, the first term of (21) becomes Finally, we obtain 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 × ( Ĥ × n ) ⋅ ( φ 0 where 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.

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 (where ∫   specifies that one that uses a Gauss quadrature rule to compute the integral) which is conserved during the time, that is, E +1 ℎ = E  ℎ , 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 E  ℎ is conserved during the discrete time: 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] Now, we examine a face Γ =  + ∩  − .We obtain the terms We can decompose  under the form 1 ± + 2 + 3, where 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].

R𝐾 ((𝑙, 𝐼) , (𝑙
The condition given by for all  is sufficient to ensure the stability of the DG scheme; where (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].

Numerical Results.
In this section, we present three examples to quantify the advantages of the proposed mixedspatial 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.he 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  7 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.
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 (y) cos () , where  =  0 √ () 2 + () 2 ,  = 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 mixedspatial 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 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.
In the last example we propose to evaluate the field scattered by an aircraft illuminated with a plane wave, by using 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.
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).

Conclusion
In this paper we presented a p-strategy with a local timestepping 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.

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

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

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

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

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

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

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

Figure 15 :
Figure 15: Test-point evaluated on the aircraft.

Figure 16 :Figure 17 :
Figure 16: Comparison on test point  for the electric field.
we consider   as the unit outward normal to the cell , we have∀ ∈ Γ,   ∘   ( x) = Jump terms associated to , by considering  = 3 and   = 2.

Table 2 :
Results obtained for the aircraft example.