Journal of Computational Methods in Physics

Volume 2013 (2013), Article ID 563480, 13 pages

http://dx.doi.org/10.1155/2013/563480

Research Article

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

Benoit Mallet,¹ Xavier Ferrieres,¹ Sebastien Pernet,¹ Jean-Baptiste Laurent,¹ Bernard Pecqueux,² and Pierre Seimandi²

¹ONERA, The French Aerospace Lab, 31055 Toulouse, France

²CEA DAM, GRAMAT, 46500 Gramat, France

Received 29 March 2013; Revised 20 June 2013; Accepted 20 June 2013

Academic Editor: Ivan D. Rukhlenko

Copyright © 2013 Benoit Mallet et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 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 of hexahedral elements being a partition of . We introduce the following space of approximation:

where 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

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 ,(ii)for all .

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 , 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 , be a set of points of , where , and are Gauss quadrature points on . At the point , we define on three basis functions , where is the Lagrange interpolation polynomial and denotes the classical cartesian base. On , the corresponding basis functions are defined by , where . Finally, the set of basis function of is

and the dimension of is (where is the number of cells). So, each can be written in this way: for all

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 ,   and ,

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 ):

(ii)for stiffness matrices, ( denotes the component of vector ):

(iii)for jump matrices

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:

where , and are 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 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 , with 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 where 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 . In other terms, the length 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 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 , for a approximation we obtain 0.1 wavelength for the ratio of the size of the cells, for we obtain 0.2, and for 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 between each degree of freedom. The mixed-spatial order strategy proposed is the following:

(i)define the smallest wavelength by considering the spectrum of the source; (ii)define an equivalent geometric length on each cell; (iii)allocate a spatial order on each cell by considering , , 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 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 and a dissipative term , where . For the evaluation of , if we consider as the unit outward normal to the cell , we have

where . The term can be written as

where

This term is also equal to

Then we obtain

Concerning the term , 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

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 , 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 on and a similar development as for the jump terms of the cell , we obtain

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.

Figure 3: Jump terms by considering .

Figure 4: Jump terms associated to , by considering and .

Figure 5: Jump terms associated to , by considering and .

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

(where specifies that one that uses a Gauss quadrature rule to compute the integral) which is conserved during the time, that is, , 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:

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 and 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 , where

So, if one takes for the Gauss quadrature rule of order and for the same quadrature rule for and (not necessarily the one corresponding to the highest spatial order), one obtains .

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

Theorem 2. Let be the matrices defined by the following: for all and for all ,

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

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

where , , and .

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 and cells projected on a plane.

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

Table 1: Time CPU and memory storage comparison between , , and mixed-spatial approaches.

Order

Mixed

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 · 
7.74  · 
5.17 · 

(a)

(b)

Figure 13: Comparison between , , 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 , we have 87.5% of cells and 12.5% of 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 in order to obtain points per wavelength for the spatial discretization. So, theoretically, one must use a 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 spatial approximation could be sufficient in some parts of the mesh. Then by using a mixed 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 , , and DG approximation for the same mesh. We can see that the solutions obtained by the mixed-spatial order is similar to the DG approximation, with a best cost in CPU-time and memory storage (see Table 2).

Table 2: Results obtained for the aircraft example.

Order

-

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 · 
2.18 · 
1.53 · 

Maximal time step
8.51 · 
5.88 · 
4.15 · 

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

References

G. Cohen, X. Ferrieres, and S. Pernet, “A spatial high-order hexahedral discontinuous Galerkin method to solve Maxwell's equations in time domain,” Journal of Computational Physics, vol. 217, no. 2, pp. 340–363, 2006.

E. Montseny, S. Pernet, X. Ferriéres, and G. Cohen, “Dissipative terms and local time-stepping improvements in a spatial high order Discontinuous Galerkin scheme for the time-domain Maxwell's equations,” Journal of Computational Physics, vol. 227, no. 14, pp. 6795–6820, 2008.

M. Dumbser, M. Käser, and E. F. Toro, “An arbitrary high-order Discontinuous Galerkin method for elastic waves on unstructured meshes—V. Local time stepping and p-adaptivity,” Geophysical Journal International, vol. 171, no. 2, pp. 695–717, 2007.

A. Taube, M. Dumbser, C.-D. Munz, and R. Schneider, “A high-order discontinuous Galerkin method with time-accurate local time stepping for the Maxwell equations,” International Journal of Numerical Modelling, vol. 22, no. 1, pp. 77–103, 2009.

S. Pernet, Étude de méthodes d'ordre élevé pour résoudre les équations de Maxwell dans le domaine temporel. Application à la détection et à la compatibilité électromagnétique [Ph.D. thesis], Université de Paris IX, 2004.

J.-P. Berenger, “Three-dimensional perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics, vol. 127, no. 2, pp. 363–379, 1996.

G. Cohen and P. Monk, “Mur-Nédélec finite element schemes for Maxwell's equations,” Computer Methods in Applied Mechanics and Engineering, vol. 169, no. 3-4, pp. 197–217, 1999.

G. C. Cohen, Higher-Order Numerical Methods for Transient Wave Equations, Springer, Berlin, Germany, 2002.

K. S. Yee, “Numerical solution of initial boundary value problem involving Maxwell's equations in isotropic media,” IEEE Transactions on Antennas and Propagation, vol. AP-14, no. 3, pp. 302–307, 1966.

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House, 3rd edition, 2005.