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

Authors: Benoit Mallet; Xavier Ferrieres; Sebastien Pernet; Jean-Baptiste Laurent; Bernard Pecqueux; Pierre Seimandi
Journal: Journal of Computational Methods in Physics
 (2013)

---

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


---

## Body

## 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 ℝ3 whose boundary is ∂Ω, and let n denote the unit outward normal to Ω. Let ɛ(x), μ(x), and σ(x) denote, respectively, the permittivity, the permeability, and the conductivity of the medium.We consider the problem described by the Maxwell equations.find(E,H):Ω×]0,T[→ℝ3×ℝ3 such that
(1)ɛ∂tE+σE=∇×H,μ∂tH=-∇×E,n×E=0on∂Ω,
where E and H 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 (Ki)i=1⋯N being a partition of Ω. We introduce the following space of approximation:
(2)Vr={v∈[L2(Ω)]3;∀K∈𝒯,DFK*v|K∘FK∈[Qr(K^)]3},
where K^=[0,1]3 is the unit cube or the reference element, for all K∈𝒯,FK:K^→K denotes the trilinear mapping which associates the vertices of each element, Qr(K^) is the space of polynomials of degree at most equal to r∈ℕ* in each variable on K^, and DFK and JK are, respectively, the Jacobian matrix and its determinant associated with the map FK. Moreover, to each K∈𝒯, we associate the outward unit normal nK.Usually, for DG methods,E and H fields are approximated by polynomials on each cell. In our case, we approximate by polynomials the fields DFK*E∘FK and DFK*H∘FK 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(Eh(·,t),Hh(·,t))∈Vr×Vr such that for all K∈𝒯 and for all ψ,ϕ∈Vr(3)∫K‍ɛ∂tEh·ψdx+∫K‍σEh·ψdx=∫K‍∇×Hh·ψdx+∫∂K‍(α∂KK⟦nK×(Eh×nK)⟧∂KK+β∂KK⟦Hh×nK⟧∂KK)·ψds,∫K‍μ∂tHh·ϕdx=-∫K‍∇×Eh·ϕdx+∫∂K‍(γ∂KK⟦Eh×nK⟧∂KK+δ∂KK⟦nK×(Hh×nK)⟧∂KK)·ϕds,
where α∂KK, β∂KK, γ∂KK, δ∂KK are parameters constant per face and ⟦v⟧ΓK=(v|K′)Γ-(v|K)Γ is the jump across the boundary Γ=K′∩K. When Γ is a boundary face (i.e., Γ=K∩∂Ω) then K′ does not exist and we define ⟦v⟧ΓK=-(v|K)Γ.The coefficientsα∂KK,β∂KK,γ∂KK,δ∂KK 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Γ=K∩K′,αΓK=αΓK′=0,δΓK=δΓK′=0,βΓK=βΓK′=-1/2,γΓK=γΓK′=1/2,(ii)
for allΓ=K∩∂Ω,βΓK=0=δΓK=αΓK,γΓK=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 nΔt and the magnetic fields at the time (n+1/2)Δt, with Δt defining the time step and n the current iteration in time.In order to define a setℬ of basis functions of Vr, we first define a set of basis functions of K^. Let x^ijk=(x^i,y^j,z^k), 1⩽i,j,k⩽r+1 be a set of points of K^, where x^i,y^j, and z^k are Gauss quadrature points on [0,1]. At the point x^ijk, we define on K^ three basis functions φ^ijkl(x^,y^,z^)=Li(x^)Lj(y^)Lk(z^)el, where Li(x^)=∏m=1,m≠ir+1‍(x^-x^m)/(∏m=1,m≠ir+1‍(x^i-x^m)) is the Lagrange interpolation polynomial and (el)l=1,2,3 denotes the classical cartesian base. On K∈𝒯, the corresponding basis functions are defined by φijkl,K∘FK(x^)=(DFK*)-1φ^ijkl(x^), where x^=(x^,y^,z^). Finally, the set of basis function of Vr is
(4)ℬ={φijkl,K:K∈𝒯,l=1,2,3,i,j,k=1,…,r+1}
and the dimension of Vr is 3(r+1)3N (where N is the number of cells). So, each U(·,t)∈Vr can be written in this way: for all K∈𝒯(5)U|K(x,t)=∑l=13‍∑ijk=1r+1‍U(t)ijkl,Kφijkl,K(x),
where U(t)ijkl,K∈ℝ is the degree of freedom associated to the basis function φijkl,K.In the sequel, we say the DG method is aQr approximation when we choose a spatial approximation order of r.Thanks to the chosen approximation space, we have for allK∈𝒯,  l,s=1,2,3 and i,j,k,m,n,p=1,…,r+1,
(6)∫K‍(∇×φijkl,K)·φmnps,Kdx=sign(JK)∫K^‍(∇^×φ^ijkl)·φ^mnpsdx^,∫∂K‍(φijkl,K×nK)·φmnps,Kds=sign(JK)∫∂K^‍(φ^ijkl×n^)·φ^mnpsds^,
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, ((M)pl denotes the (p,l) component of the matrix M):
(7)∫K‍ɛU·φijkl,Kdx=ωijk∑p=13‍Uijkp,K|JK|(DFK-1(ɛ∘FK)DFK*-1)pl(x^ijk),(ii)
for stiffness matrices, ((M)l denotes the component l of vector M):
(8)∫K‍∇×U·φijkl,K=sign(JK)∑mnq=1r+1‍∑p=13‍ωijkUijkl,K(∇^×φ^ijkm(x^ijk))l,(iii)
for jump matrices(9)∫∂K‍⟦H×nK⟧∂K·φijkl,Kds=-sign(JK′)∫∂K^‍(DFK-1∘SK′KDFK′)×(H^K′×n^)·φ^ijkl∘SK′Kds^+sign(JK)∫∂K^‍(H^K×n^)·φ^ijklds^.In the above expressions ωijk is the quadrature weight at point x^ijk and SK′K=FK-1∘FK′. DFK-1∘SK′KDFK′ is a permutation matrix constant per face [1].The DG formulation (3) finally leads to the semidiscrete formulation:
(10)Mɛ∂tE+MσE=RH-SiH,Mμ∂tH=-RE+SiE+SbE,
where Mɛ,Mμ, and Mσ are 3×3 block-diagonal matrices, R, the stiffness matrix, and Si and Sb 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 K. Stiffness and jump matrices just require to store the sign of the Jacobians JK 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.


### 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 N classes of cells where each class i has a time step given by 3N-iΔtminr, with Δtmin 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 Figure1. 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.


## 2.1. Discontinuous Galerkin Method

LetΩ be a bounded open subset of ℝ3 whose boundary is ∂Ω, and let n denote the unit outward normal to Ω. Let ɛ(x), μ(x), and σ(x) denote, respectively, the permittivity, the permeability, and the conductivity of the medium.We consider the problem described by the Maxwell equations.find(E,H):Ω×]0,T[→ℝ3×ℝ3 such that
(1)ɛ∂tE+σE=∇×H,μ∂tH=-∇×E,n×E=0on∂Ω,
where E and H 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 (Ki)i=1⋯N being a partition of Ω. We introduce the following space of approximation:
(2)Vr={v∈[L2(Ω)]3;∀K∈𝒯,DFK*v|K∘FK∈[Qr(K^)]3},
where K^=[0,1]3 is the unit cube or the reference element, for all K∈𝒯,FK:K^→K denotes the trilinear mapping which associates the vertices of each element, Qr(K^) is the space of polynomials of degree at most equal to r∈ℕ* in each variable on K^, and DFK and JK are, respectively, the Jacobian matrix and its determinant associated with the map FK. Moreover, to each K∈𝒯, we associate the outward unit normal nK.Usually, for DG methods,E and H fields are approximated by polynomials on each cell. In our case, we approximate by polynomials the fields DFK*E∘FK and DFK*H∘FK 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(Eh(·,t),Hh(·,t))∈Vr×Vr such that for all K∈𝒯 and for all ψ,ϕ∈Vr(3)∫K‍ɛ∂tEh·ψdx+∫K‍σEh·ψdx=∫K‍∇×Hh·ψdx+∫∂K‍(α∂KK⟦nK×(Eh×nK)⟧∂KK+β∂KK⟦Hh×nK⟧∂KK)·ψds,∫K‍μ∂tHh·ϕdx=-∫K‍∇×Eh·ϕdx+∫∂K‍(γ∂KK⟦Eh×nK⟧∂KK+δ∂KK⟦nK×(Hh×nK)⟧∂KK)·ϕds,
where α∂KK, β∂KK, γ∂KK, δ∂KK are parameters constant per face and ⟦v⟧ΓK=(v|K′)Γ-(v|K)Γ is the jump across the boundary Γ=K′∩K. When Γ is a boundary face (i.e., Γ=K∩∂Ω) then K′ does not exist and we define ⟦v⟧ΓK=-(v|K)Γ.The coefficientsα∂KK,β∂KK,γ∂KK,δ∂KK 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Γ=K∩K′,αΓK=αΓK′=0,δΓK=δΓK′=0,βΓK=βΓK′=-1/2,γΓK=γΓK′=1/2,(ii)
for allΓ=K∩∂Ω,βΓK=0=δΓK=αΓK,γΓK=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 nΔt and the magnetic fields at the time (n+1/2)Δt, with Δt defining the time step and n the current iteration in time.In order to define a setℬ of basis functions of Vr, we first define a set of basis functions of K^. Let x^ijk=(x^i,y^j,z^k), 1⩽i,j,k⩽r+1 be a set of points of K^, where x^i,y^j, and z^k are Gauss quadrature points on [0,1]. At the point x^ijk, we define on K^ three basis functions φ^ijkl(x^,y^,z^)=Li(x^)Lj(y^)Lk(z^)el, where Li(x^)=∏m=1,m≠ir+1‍(x^-x^m)/(∏m=1,m≠ir+1‍(x^i-x^m)) is the Lagrange interpolation polynomial and (el)l=1,2,3 denotes the classical cartesian base. On K∈𝒯, the corresponding basis functions are defined by φijkl,K∘FK(x^)=(DFK*)-1φ^ijkl(x^), where x^=(x^,y^,z^). Finally, the set of basis function of Vr is
(4)ℬ={φijkl,K:K∈𝒯,l=1,2,3,i,j,k=1,…,r+1}
and the dimension of Vr is 3(r+1)3N (where N is the number of cells). So, each U(·,t)∈Vr can be written in this way: for all K∈𝒯(5)U|K(x,t)=∑l=13‍∑ijk=1r+1‍U(t)ijkl,Kφijkl,K(x),
where U(t)ijkl,K∈ℝ is the degree of freedom associated to the basis function φijkl,K.In the sequel, we say the DG method is aQr approximation when we choose a spatial approximation order of r.Thanks to the chosen approximation space, we have for allK∈𝒯,  l,s=1,2,3 and i,j,k,m,n,p=1,…,r+1,
(6)∫K‍(∇×φijkl,K)·φmnps,Kdx=sign(JK)∫K^‍(∇^×φ^ijkl)·φ^mnpsdx^,∫∂K‍(φijkl,K×nK)·φmnps,Kds=sign(JK)∫∂K^‍(φ^ijkl×n^)·φ^mnpsds^,
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, ((M)pl denotes the (p,l) component of the matrix M):
(7)∫K‍ɛU·φijkl,Kdx=ωijk∑p=13‍Uijkp,K|JK|(DFK-1(ɛ∘FK)DFK*-1)pl(x^ijk),(ii)
for stiffness matrices, ((M)l denotes the component l of vector M):
(8)∫K‍∇×U·φijkl,K=sign(JK)∑mnq=1r+1‍∑p=13‍ωijkUijkl,K(∇^×φ^ijkm(x^ijk))l,(iii)
for jump matrices(9)∫∂K‍⟦H×nK⟧∂K·φijkl,Kds=-sign(JK′)∫∂K^‍(DFK-1∘SK′KDFK′)×(H^K′×n^)·φ^ijkl∘SK′Kds^+sign(JK)∫∂K^‍(H^K×n^)·φ^ijklds^.In the above expressions ωijk is the quadrature weight at point x^ijk and SK′K=FK-1∘FK′. DFK-1∘SK′KDFK′ is a permutation matrix constant per face [1].The DG formulation (3) finally leads to the semidiscrete formulation:
(10)Mɛ∂tE+MσE=RH-SiH,Mμ∂tH=-RE+SiE+SbE,
where Mɛ,Mμ, and Mσ are 3×3 block-diagonal matrices, R, the stiffness matrix, and Si and Sb 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 K. Stiffness and jump matrices just require to store the sign of the Jacobians JK 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.


## 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 N classes of cells where each class i has a time step given by 3N-iΔtminr, with Δtmin 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 Figure1. 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/n where λmin defines the smallest wavelength of the source and n is an integer which depends on the wanted accuracy. For example, in an FDTD approach, generally, people use n=10. In other terms, the length λmin/n 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 r, to keep a discretization of λmin/n between each degree of freedom, the ratio of the size of the cell by the smallest wavelength will be given by r/n. For example, by considering n=10, for a Q1 approximation we obtain 0.1 wavelength for the ratio of the size of the cells, for Q2 we obtain 0.2, and for Q3 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/n 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 lengthh on each cell;(iii)
allocate a spatial orderr=(h/λmin)n on each cell by considering h, λmin, and a given n.The difficulty in this strategy is to determine a value ofh 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 Section2. 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 φi0j0k0l0,K be a given basis function on K; by considering two contiguous cells K and K′ with different spatial approximated orders r and r′, r>r′, the jump term for the equation of the magnetic field in the system (3) on the cell K is split into a centered term ScK=∫Γ‍⟦E×nK⟧·φi0j0k0l0,Kds and a dissipative term SdK=∫Γ‍⟦nK×(H×nK)⟧·φi0j0k0l0,Kds, where Γ=K∩K′. For the evaluation of ScK=∫Γ‍(EK′×nK)·φi0,j0,k0l0,Kds-∫Γ‍(EK×nK)·φi0j0k0l0,Kds, if we consider nK as the unit outward normal to the cell K, we have
(11)∀x∈Γ,nK∘FK(x^)=DFK*-1(x^)n^∥DFK*-1(x^)n^∥,ds=∥DFK*-1(x^)n^∥|JK(x^)|,
where x=FK(x^). The term EK(x)×nK can be written as
(12)DFK*-1(x^)E^(x^)×DFK*-1(x^)n^∥DFK*-1(x^)n^∥,
where
(13)E^(x^)=∑l=13‍∑ijk=1r+13‍Eijkl,Kφ^ijkl(x^).
This term is also equal to
(14)DFK(x^)JK(x^)E^(x^)×n^.
Then we obtain
(15)∫Γ‍(EK×nK)·φi0j0k0l0,Kds=∫Γ^‍sign(JK(x^))(E^(x^)×n^)·φ^i0j0k0l0(x^)ds^.
Concerning the term ∫Γ‍(EK′×nK)·φi0j0k0l0,Kds, we take SKK′=FK′-1∘FK with FK′ being the trilinear transformation between the cell K′ and the reference element K^. Then, by using nK=-nK′, the previous integral term can be written as
(16)-∫Γ^‍DFK′∘SKK′(x^)JK′∘SKK′(x^)E^K′×n^′∥DFK′*-1(x^)n^′∥·DFK*-1φ^i0j0k0l0(x^)∥DFK*-1(x^)n^∥|JK(x^)|ds^
with n′ defined as
(17)n^K′=DFK′*-1∘SKK′(x^)n^′∥DFK′*-1∘SKK′(x^)n^′∥.
By using the relation
(18)∥DFK*-1(x^)n^∥|JK(x^)|=∥DFK′*-1∘SKK′(x^)n^′∥|JK′∘SKK′(x^)|
we obtain for the previous term
(19)-∫Γ^‍sign(JK∘SKK′)(x^)DFK-1(x^)DFK′∘SKK′(x^)E^K′(x^)×n^′·φ^i0,j0,k0l0(x^)ds^.
Finally, for the term ScK, we have(20)ScK=-∫Γ^‍sign(JK∘SKK′)(x^)DFK-1(x^)DFK′∘SKK′(x^)E^K′(x^)×n^′·φ^i0,j0,k0l0(x^)ds^-∫Γ^‍sign(JK(x^))(E^(x^)×n^)·φ^i0j0k0l0(x^)ds^.For the dissipative term we have(21)SdK=∫Γ‍nK×(HK′×nK)·φi0j0k0l0,Kds-∫Γ‍nK×(HK×nK)·φi0j0k0l0,Kds.
By using H^(x^)=∑l=13‍∑ijk=1(r+1)3Hijkl,Kφ^ijkl(x^), the second term of (21) can be written as
(22)∫Γ‍(HK×nK)·φi0j0k0l0,K×nKds=∫Γ^‍DFK(x^)JK(x^)H^K(x^)×n^∥DFK*-1(x^)n^∥·DFK(x^)JK(x^)(φ^i0j0k0l0×n^)∥DFK*-1(x^)n^∥|JK(x^)|ds^=∫Γ^‍sign(JK(x^))DFK*(x^)DFK(x^)JK(x^)∥DFK*-1(x^)n^∥×(H^K×n^)·(φ^i0j0k0l0×n^)ds^.
In the same way, the first term of (21) becomes
(23)∫Γ^‍DFK(x^)DFK′∘SKK′(x^)JK′∘SKK′(x^)∥DFK′*-1∘SKK′(x^)n^′∥×(H^K′×n^′)·(φ^i0j0k0l0×n^)ds^.
Finally, we obtain
(24)SdK=∫Γ^‍DFK(x^)DFK′∘SKK′(x^)JK′∘SKK′(x^)∥DFK′*-1∘SKK′(x^)n^′∥×(H^K′×n^′)·(φ^i0j0k0l0×n^)ds^-∫Γ^‍sign(JK(x^))DFK*(x^)DFK(x^)JK(x^)∥DFK*-1(x^)n^∥×(H^K×n^)·(φ^i0j0k0l0×n^)ds^.Similarly, the jump term on the cellK′ is given by a centered term ScK′=∫Γ‍⟦E×n′⟧·φ′ds and a dissipative term SdK′=∫Γ‍⟦n′×(H×n′)⟧·φ′ds. By considering a basis function φi0,j0,k0l0,K′ on K′ and a similar development as for the jump terms of the cell K, we obtain
(25)ScK′=-∫Γ^‍sign(JK′∘SK′K)(x^)DFK′-1(x^)DFK∘SK′K(x^)E^K(x^)×n^′·φ^i0,j0,k0l0(x^)ds^-∫Γ^‍sign(JK′(x^))(E^K′(x^)×n^)·φ^i0j0k0l0(x^)ds^,SdK′=∫Γ^‍DFK′(x^)DFK∘SK′K(x^)JK∘SK′K(x^)∥DFK*-1∘SK′K(x^)n^′∥×(H^K×n^′)·(φ^i0j0k0l0×n^)ds^-∫Γ^‍sign(JK′(x^))DFK′*(x^)DFK′(x^)JK′(x^)∥DFK′*-1(x^)n^∥×(H^K′×n^)·(φ^i0j0k0l0×n^)ds^,
where
(26)Γ=FK′(Γ^)=FK(Γ^′),SK′K=FK*-1∘FK′,x^=FK′-1(x),nK′∘FK′(x^)=DFK′*-1(x^)n^∥DFK′*-1(x^)n^∥,nK∘FK′(x^)=DFK*-1∘SK′K(x^)n^′∥DFK*-1∘SK′K(x^)n^′∥.By using a Gauss quadrature formula of order equal to the largest spatial approximation order ofK and K′ 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 consideringr=r′=3.Figure 4
Jump terms associated toK, by considering r=3 and r′=2.Figure 5
Jump terms associated toK′, by considering r=3 and r′=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(27)ℰhn=∑K∈𝒯h‍(∫KG‍ɛKEh,Kn·Eh,Kndx+∫KG‍μKHh,Kn+1/2·Hh,Kn-1/2dx)
(where ∫KG‍ specifies that one that uses a Gauss quadrature rule to compute the integral) which is conserved during the time, that is, ℰhn+1=ℰhn, 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ℰhn is conserved during the discrete time:
(28)∀n≥0,ℰhn+1-ℰhn=0.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ψ=Ehn+1+Ehn and ϕ=Hhn+1/2 in the nondissipative fully discrete scheme, one can easily obtain [5]
(29)ℰhn+1-ℰhn=∑K∈𝒯‍[∫∂KG‍γ⟦(Ehn+1+Ehn)×nK⟧∂KK∫K‍∇×HKn+1/2·(EKn+1+EKn)dx∑K∈𝒯‍-∫K‍∇×(EKn+1+EKn)·HKn+1/2dx∑K∈𝒯‍+∫∂KG‍β⟦Hhn+1/2×nK⟧∂KK·(EKn+1+EKn)∑K∈𝒯‍+∫∂KG‍γ⟦(Ehn+1+Ehn)×nK⟧∂KK·HKn+1/2]=∑K∈𝒯‍[∫∂K‍(nK×HKn+1/2)·(EKn+1+EKn)dσ∑K∈𝒯‍+∫∂KG‍β⟦Hhn+1/2×nK⟧∂KK·(EKn+1+EKn)∑K∈𝒯‍+∫∂KG‍γ⟦(Ehn+1+Ehn)×nK⟧∂KK·HKn+1/2].
Now, we examine a face Γ=K+¯∩K-¯. We obtain the terms
(30)T=∫Γ‍[n+×H+n+1/2·(E+n+1+E+n)+n-×H-n+1/2·(E-n+1+E-n)]dσ+∫ΓG+‍β(H-n+1/2×n+-H+n+1/2×n+)·(E+n+1+E+n)+∫ΓG-‍β(H+n+1/2×n--H-n+1/2×n-)·(E-n+1+E-n)+∫ΓG+‍γ((E-n+1+E-n)×n+-(E+n+1+E+n)×n+)·H+n+1/2+∫ΓG-‍γ((E+n+1+E+n)×n--(E-n+1+E-n)×n-)·H-n+1/2.
We can decompose T under the form T1±+T2+T3, where
(31)T1±=∫Γ‍n±×H±n+1/2·(E±n+1+E±n)dσ-∫ΓG±‍βH±n+1/2×n±·(E±n+1+E±n)-∫ΓG±‍γ(E±n+1+E±n)×n±·H±n+1/2,T2=∫ΓG+‍βH-n+1/2×n+·(E+n+1+E+n)-∫ΓG-‍γ(E+n+1+E+n)×n-·H-n+1/2,T3=∫ΓG-‍βH+n+1/2×n-·(E-n+1+E-n)-∫ΓG+‍γ(E-n+1+E-n)×n+·H+n+1/2.
So, if one takes forT1± the Gauss quadrature rule of order r± and for T2,T3 the same quadrature rule for G+ and G- (not necessarily the one corresponding to the highest spatial order), one obtains T=0.So, one can give the stability result [1].Theorem 2.
Letℛ^K,𝒟^K,ℬ^K be the 3(rK+1)3×3(rK+1)3 matrices defined by the following: for all l,l′∈{1,2,3} and for all I=(i,j,k), I′=(i′,j′,k′)∈{1,…,rK+1}3(32)ℛ^K((l,I),(l′,I′))=∫∂K^‍∇^×φ^ijkl·∇^×φ^i′j′k′l′,𝒟^K((l,I),(l′,I′))=δll′δII′ωijk,ℬ^K((l,I),(l′,I′))=∫∂K^‍(φ^ijkl×n^)·(φ^i′j′k′l′×n^)ds^.
The condition given by
(33)ΔtΛK<2cK1×((12max1≤i≤nbfiK(μKμV(i,K),ɛKɛV(i,K)))λmax(𝒟^K-1/2ℛ^K𝒟^K-1/2)+(12max1≤i≤nbfiK(μKμV(i,K),ɛKɛV(i,K)))×λmax(𝒟^K-1/2ℬ^K𝒟^K-1/2)(12max1≤i≤nbfiK(μKμV(i,K),ɛKɛV(i,K))))-1
for all K is sufficient to ensure the stability of the DG scheme; where
(34)ΛK=min1≤i,j,k≤rK+1(|JK(x^ijk)|λmax((DFK*DFK)(x^ijk)))
(in the case of a uniform cartesian grid whose spatial step is h, ΛK=h), where nbfiK is the number of internal faces of K and V(i,K) is the the neighbor of K.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 a1m×1m 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 r=7. So, we obtain a Q7 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: a3×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
(35)Ex=0,Ey=0,Ez=sin(mπx)sin(nπy)cos(ωt),Hx=-1ωμ0πnsin(mπx)cos(nπy)sin(ωt),Hy=1ωμ0πmcos(mπx)sin(nπy)sin(ωt),Hz=0,
where ω=c0(mπ)2+(nπ)2, m=1, and n=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 Figure12.Figure 12
Distributions forQ3 and Q4 cells projected on a (x,y) plane.Figure13 shows the analytical solution and the solutions obtained by the Q3, Q4, and the mixed Q3-Q4 DG approximations. We can see that the Q3-Q4 DG approximation leads to the same accurate solution of the Q4 DG method but its numerical costs are far from smaller (see Table 1). Finally, we show in Figure 14 that the Q3 solution is not sufficiently accurate whereas the mixed Q3-Q4 approximation gives the best compromise between CPU-time, memory storage, and accuracy.Table 1
Time CPU and memory storage comparison betweenQ3, Q4, and mixed-spatial Q3-Q4 approaches.
Order
Q
3
MixedQ3-Q4
Q
4
Ratio
100%
87.5%–12.5%
100%
CPU time (s)
25009
44217
70710
Memory storage (MB)
23
25
35
DOF
294912
330048
576000
Time step
7.74 ·10-12
7.74  ·10-12
5.17 ·10-12Comparison betweenQ3, Q4, mixed-spatial order, and analytic solutions for the electric field Ez 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 Q3-Q4, we have 87.5% of Q3 cells and 12.5% of Q4 cells.
(a)
(b)Point to point error for the different computed solutions and L2 error.
(a)
(b)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 Section2.2, jointly with the mixed-spatial order strategy.By considering the mesh (see Figure15) 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 Q2 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 Q1 spatial approximation could be sufficient in some parts of the mesh. Then by using a mixed Q1-Q2 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 A (see Figure 15) obtained by Q1, Q2, and Q1-Q2 DG approximation for the same mesh. We can see that the solutions obtained by the mixed-spatial order is similar to the Q2 DG approximation, with a best cost in CPU-time and memory storage (see Table 2).Table 2
Results obtained for the aircraft example.
Order
Q
1
Q
1-Q2
Q
2
Ratio
100%
20%–80%
100%
CPU time (s)
16478
58445
76229
Memory storage (MB)
673
1208
1316
DOF
5150208
14918982
17381952
Minimal time step
3.15 ·10-12
2.18 ·10-12
1.53 ·10-12
Maximal time step
8.51 ·10-11
5.88 ·10-11
4.15 ·10-11Figure 15
Test-point evaluated on the aircraft.Comparison on test pointA for the electric field.
(a)
(b)Comparison on test pointA for the magnetic field.
(a)
(b)


## 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 Section2. 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 φi0j0k0l0,K be a given basis function on K; by considering two contiguous cells K and K′ with different spatial approximated orders r and r′, r>r′, the jump term for the equation of the magnetic field in the system (3) on the cell K is split into a centered term ScK=∫Γ‍⟦E×nK⟧·φi0j0k0l0,Kds and a dissipative term SdK=∫Γ‍⟦nK×(H×nK)⟧·φi0j0k0l0,Kds, where Γ=K∩K′. For the evaluation of ScK=∫Γ‍(EK′×nK)·φi0,j0,k0l0,Kds-∫Γ‍(EK×nK)·φi0j0k0l0,Kds, if we consider nK as the unit outward normal to the cell K, we have
(11)∀x∈Γ,nK∘FK(x^)=DFK*-1(x^)n^∥DFK*-1(x^)n^∥,ds=∥DFK*-1(x^)n^∥|JK(x^)|,
where x=FK(x^). The term EK(x)×nK can be written as
(12)DFK*-1(x^)E^(x^)×DFK*-1(x^)n^∥DFK*-1(x^)n^∥,
where
(13)E^(x^)=∑l=13‍∑ijk=1r+13‍Eijkl,Kφ^ijkl(x^).
This term is also equal to
(14)DFK(x^)JK(x^)E^(x^)×n^.
Then we obtain
(15)∫Γ‍(EK×nK)·φi0j0k0l0,Kds=∫Γ^‍sign(JK(x^))(E^(x^)×n^)·φ^i0j0k0l0(x^)ds^.
Concerning the term ∫Γ‍(EK′×nK)·φi0j0k0l0,Kds, we take SKK′=FK′-1∘FK with FK′ being the trilinear transformation between the cell K′ and the reference element K^. Then, by using nK=-nK′, the previous integral term can be written as
(16)-∫Γ^‍DFK′∘SKK′(x^)JK′∘SKK′(x^)E^K′×n^′∥DFK′*-1(x^)n^′∥·DFK*-1φ^i0j0k0l0(x^)∥DFK*-1(x^)n^∥|JK(x^)|ds^
with n′ defined as
(17)n^K′=DFK′*-1∘SKK′(x^)n^′∥DFK′*-1∘SKK′(x^)n^′∥.
By using the relation
(18)∥DFK*-1(x^)n^∥|JK(x^)|=∥DFK′*-1∘SKK′(x^)n^′∥|JK′∘SKK′(x^)|
we obtain for the previous term
(19)-∫Γ^‍sign(JK∘SKK′)(x^)DFK-1(x^)DFK′∘SKK′(x^)E^K′(x^)×n^′·φ^i0,j0,k0l0(x^)ds^.
Finally, for the term ScK, we have(20)ScK=-∫Γ^‍sign(JK∘SKK′)(x^)DFK-1(x^)DFK′∘SKK′(x^)E^K′(x^)×n^′·φ^i0,j0,k0l0(x^)ds^-∫Γ^‍sign(JK(x^))(E^(x^)×n^)·φ^i0j0k0l0(x^)ds^.For the dissipative term we have(21)SdK=∫Γ‍nK×(HK′×nK)·φi0j0k0l0,Kds-∫Γ‍nK×(HK×nK)·φi0j0k0l0,Kds.
By using H^(x^)=∑l=13‍∑ijk=1(r+1)3Hijkl,Kφ^ijkl(x^), the second term of (21) can be written as
(22)∫Γ‍(HK×nK)·φi0j0k0l0,K×nKds=∫Γ^‍DFK(x^)JK(x^)H^K(x^)×n^∥DFK*-1(x^)n^∥·DFK(x^)JK(x^)(φ^i0j0k0l0×n^)∥DFK*-1(x^)n^∥|JK(x^)|ds^=∫Γ^‍sign(JK(x^))DFK*(x^)DFK(x^)JK(x^)∥DFK*-1(x^)n^∥×(H^K×n^)·(φ^i0j0k0l0×n^)ds^.
In the same way, the first term of (21) becomes
(23)∫Γ^‍DFK(x^)DFK′∘SKK′(x^)JK′∘SKK′(x^)∥DFK′*-1∘SKK′(x^)n^′∥×(H^K′×n^′)·(φ^i0j0k0l0×n^)ds^.
Finally, we obtain
(24)SdK=∫Γ^‍DFK(x^)DFK′∘SKK′(x^)JK′∘SKK′(x^)∥DFK′*-1∘SKK′(x^)n^′∥×(H^K′×n^′)·(φ^i0j0k0l0×n^)ds^-∫Γ^‍sign(JK(x^))DFK*(x^)DFK(x^)JK(x^)∥DFK*-1(x^)n^∥×(H^K×n^)·(φ^i0j0k0l0×n^)ds^.Similarly, the jump term on the cellK′ is given by a centered term ScK′=∫Γ‍⟦E×n′⟧·φ′ds and a dissipative term SdK′=∫Γ‍⟦n′×(H×n′)⟧·φ′ds. By considering a basis function φi0,j0,k0l0,K′ on K′ and a similar development as for the jump terms of the cell K, we obtain
(25)ScK′=-∫Γ^‍sign(JK′∘SK′K)(x^)DFK′-1(x^)DFK∘SK′K(x^)E^K(x^)×n^′·φ^i0,j0,k0l0(x^)ds^-∫Γ^‍sign(JK′(x^))(E^K′(x^)×n^)·φ^i0j0k0l0(x^)ds^,SdK′=∫Γ^‍DFK′(x^)DFK∘SK′K(x^)JK∘SK′K(x^)∥DFK*-1∘SK′K(x^)n^′∥×(H^K×n^′)·(φ^i0j0k0l0×n^)ds^-∫Γ^‍sign(JK′(x^))DFK′*(x^)DFK′(x^)JK′(x^)∥DFK′*-1(x^)n^∥×(H^K′×n^)·(φ^i0j0k0l0×n^)ds^,
where
(26)Γ=FK′(Γ^)=FK(Γ^′),SK′K=FK*-1∘FK′,x^=FK′-1(x),nK′∘FK′(x^)=DFK′*-1(x^)n^∥DFK′*-1(x^)n^∥,nK∘FK′(x^)=DFK*-1∘SK′K(x^)n^′∥DFK*-1∘SK′K(x^)n^′∥.By using a Gauss quadrature formula of order equal to the largest spatial approximation order ofK and K′ 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 consideringr=r′=3.Figure 4
Jump terms associated toK, by considering r=3 and r′=2.Figure 5
Jump terms associated toK′, by considering r=3 and r′=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(27)ℰhn=∑K∈𝒯h‍(∫KG‍ɛKEh,Kn·Eh,Kndx+∫KG‍μKHh,Kn+1/2·Hh,Kn-1/2dx)
(where ∫KG‍ specifies that one that uses a Gauss quadrature rule to compute the integral) which is conserved during the time, that is, ℰhn+1=ℰhn, 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ℰhn is conserved during the discrete time:
(28)∀n≥0,ℰhn+1-ℰhn=0.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ψ=Ehn+1+Ehn and ϕ=Hhn+1/2 in the nondissipative fully discrete scheme, one can easily obtain [5]
(29)ℰhn+1-ℰhn=∑K∈𝒯‍[∫∂KG‍γ⟦(Ehn+1+Ehn)×nK⟧∂KK∫K‍∇×HKn+1/2·(EKn+1+EKn)dx∑K∈𝒯‍-∫K‍∇×(EKn+1+EKn)·HKn+1/2dx∑K∈𝒯‍+∫∂KG‍β⟦Hhn+1/2×nK⟧∂KK·(EKn+1+EKn)∑K∈𝒯‍+∫∂KG‍γ⟦(Ehn+1+Ehn)×nK⟧∂KK·HKn+1/2]=∑K∈𝒯‍[∫∂K‍(nK×HKn+1/2)·(EKn+1+EKn)dσ∑K∈𝒯‍+∫∂KG‍β⟦Hhn+1/2×nK⟧∂KK·(EKn+1+EKn)∑K∈𝒯‍+∫∂KG‍γ⟦(Ehn+1+Ehn)×nK⟧∂KK·HKn+1/2].
Now, we examine a face Γ=K+¯∩K-¯. We obtain the terms
(30)T=∫Γ‍[n+×H+n+1/2·(E+n+1+E+n)+n-×H-n+1/2·(E-n+1+E-n)]dσ+∫ΓG+‍β(H-n+1/2×n+-H+n+1/2×n+)·(E+n+1+E+n)+∫ΓG-‍β(H+n+1/2×n--H-n+1/2×n-)·(E-n+1+E-n)+∫ΓG+‍γ((E-n+1+E-n)×n+-(E+n+1+E+n)×n+)·H+n+1/2+∫ΓG-‍γ((E+n+1+E+n)×n--(E-n+1+E-n)×n-)·H-n+1/2.
We can decompose T under the form T1±+T2+T3, where
(31)T1±=∫Γ‍n±×H±n+1/2·(E±n+1+E±n)dσ-∫ΓG±‍βH±n+1/2×n±·(E±n+1+E±n)-∫ΓG±‍γ(E±n+1+E±n)×n±·H±n+1/2,T2=∫ΓG+‍βH-n+1/2×n+·(E+n+1+E+n)-∫ΓG-‍γ(E+n+1+E+n)×n-·H-n+1/2,T3=∫ΓG-‍βH+n+1/2×n-·(E-n+1+E-n)-∫ΓG+‍γ(E-n+1+E-n)×n+·H+n+1/2.
So, if one takes forT1± the Gauss quadrature rule of order r± and for T2,T3 the same quadrature rule for G+ and G- (not necessarily the one corresponding to the highest spatial order), one obtains T=0.So, one can give the stability result [1].Theorem 2.
Letℛ^K,𝒟^K,ℬ^K be the 3(rK+1)3×3(rK+1)3 matrices defined by the following: for all l,l′∈{1,2,3} and for all I=(i,j,k), I′=(i′,j′,k′)∈{1,…,rK+1}3(32)ℛ^K((l,I),(l′,I′))=∫∂K^‍∇^×φ^ijkl·∇^×φ^i′j′k′l′,𝒟^K((l,I),(l′,I′))=δll′δII′ωijk,ℬ^K((l,I),(l′,I′))=∫∂K^‍(φ^ijkl×n^)·(φ^i′j′k′l′×n^)ds^.
The condition given by
(33)ΔtΛK<2cK1×((12max1≤i≤nbfiK(μKμV(i,K),ɛKɛV(i,K)))λmax(𝒟^K-1/2ℛ^K𝒟^K-1/2)+(12max1≤i≤nbfiK(μKμV(i,K),ɛKɛV(i,K)))×λmax(𝒟^K-1/2ℬ^K𝒟^K-1/2)(12max1≤i≤nbfiK(μKμV(i,K),ɛKɛV(i,K))))-1
for all K is sufficient to ensure the stability of the DG scheme; where
(34)ΛK=min1≤i,j,k≤rK+1(|JK(x^ijk)|λmax((DFK*DFK)(x^ijk)))
(in the case of a uniform cartesian grid whose spatial step is h, ΛK=h), where nbfiK is the number of internal faces of K and V(i,K) is the the neighbor of K.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 a1m×1m 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 r=7. So, we obtain a Q7 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: a3×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
(35)Ex=0,Ey=0,Ez=sin(mπx)sin(nπy)cos(ωt),Hx=-1ωμ0πnsin(mπx)cos(nπy)sin(ωt),Hy=1ωμ0πmcos(mπx)sin(nπy)sin(ωt),Hz=0,
where ω=c0(mπ)2+(nπ)2, m=1, and n=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 Figure12.Figure 12
Distributions forQ3 and Q4 cells projected on a (x,y) plane.Figure13 shows the analytical solution and the solutions obtained by the Q3, Q4, and the mixed Q3-Q4 DG approximations. We can see that the Q3-Q4 DG approximation leads to the same accurate solution of the Q4 DG method but its numerical costs are far from smaller (see Table 1). Finally, we show in Figure 14 that the Q3 solution is not sufficiently accurate whereas the mixed Q3-Q4 approximation gives the best compromise between CPU-time, memory storage, and accuracy.Table 1
Time CPU and memory storage comparison betweenQ3, Q4, and mixed-spatial Q3-Q4 approaches.
Order
Q
3
MixedQ3-Q4
Q
4
Ratio
100%
87.5%–12.5%
100%
CPU time (s)
25009
44217
70710
Memory storage (MB)
23
25
35
DOF
294912
330048
576000
Time step
7.74 ·10-12
7.74  ·10-12
5.17 ·10-12Comparison betweenQ3, Q4, mixed-spatial order, and analytic solutions for the electric field Ez 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 Q3-Q4, we have 87.5% of Q3 cells and 12.5% of Q4 cells.
(a)
(b)Point to point error for the different computed solutions and L2 error.
(a)
(b)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 Section2.2, jointly with the mixed-spatial order strategy.By considering the mesh (see Figure15) 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 Q2 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 Q1 spatial approximation could be sufficient in some parts of the mesh. Then by using a mixed Q1-Q2 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 A (see Figure 15) obtained by Q1, Q2, and Q1-Q2 DG approximation for the same mesh. We can see that the solutions obtained by the mixed-spatial order is similar to the Q2 DG approximation, with a best cost in CPU-time and memory storage (see Table 2).Table 2
Results obtained for the aircraft example.
Order
Q
1
Q
1-Q2
Q
2
Ratio
100%
20%–80%
100%
CPU time (s)
16478
58445
76229
Memory storage (MB)
673
1208
1316
DOF
5150208
14918982
17381952
Minimal time step
3.15 ·10-12
2.18 ·10-12
1.53 ·10-12
Maximal time step
8.51 ·10-11
5.88 ·10-11
4.15 ·10-11Figure 15
Test-point evaluated on the aircraft.Comparison on test pointA for the electric field.
(a)
(b)Comparison on test pointA for the magnetic field.
(a)
(b)


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



---

*Source: 563480-2013-07-30.xml*