Two important features of discontinuous Galerkin methods are their
flexibility with regards to the local approximation of the field
quantities and their natural ability to deal with non-conforming
meshes. The non-conformity can result from a local refinement of the
mesh (h-adaptivity), or of the approximation order (p-adaptivity)
or of both of them (hp-adaptivity). Three main reasons motivate our
investigation of DG methods able to deal with non-conforming
simplicial meshes: (1) the final mesh is made of separately
constructed meshes and a conforming assembling is a very hard task,
(2) the presence of geometrical features (fine structures, corners,
etc.) or physical ones (multi-scale problems, etc.) requiring a
highly localized refinement and, (3) the construction of hybrid meshes
mixing elements of different types such as hexahedra and tetrahedra
without the addition of another type of element (a pyramid for example
in this case). In the context of my PhD thesis, I have developed a
DGTD method based on high-order nodal basis functions for the
approximation of the electromagnetic field within a simplex, a
centered scheme for the calculation of the numerical flux at an
interface between neighboring elements, and a second-order leap-frog
time integration scheme [J1] -
[C2] -
[C3]. Moreover, the mesh is refined locally in a
non-conforming way resulting in arbitrary level hanging nodes. The
resulting method is non-dissipative, stable under some CFL-like
condition, conserves a discrete version of the electromagnetic energy,
and does not introduce much dispersion error [R1]. In this
context, a hp-like DGTD method which allows for both a local
non-conforming refinement of the mesh and a locally defined
approximation order has been designed, anlayzed and evaluated in the
context of the numerical solution of the time domain Maxwell equations
on simplicial meshes [J2] - [R3].
The space discretized discontinuous Galerkin method for the Maxwell equations
can be written in the following form
where Y=(E,H)T and H is a skew-symmetric
matrix depends only on the spatial configuration. The formal solution
of the ODE, Eq. (1), is given
where Y(0) represents the intial state of the fields and the
operator Φ(t)=etH denotes their time-evolution
matrix. Since H is skew-symmetric, the time evolution evolution
operator, Φ, is an orthogonal
transformation. Mathematically speaking, the time evolution
operator Φ(t) rotates the vector Y(t) without
changing its length ∥Y∥. In physical terms, this
means that the total electromagnetic energy is conserved, as can be
expected on physical grounds.
Seeking a time discrete solution of Eq. (1), a discretization
in time with a global time step Δt is introduced. The
construction of high-order, one-step, explicit time integrators is
based on the approximation of the matrix
exponential Φ(Δt)=eΔtH. At this
point, we invoked two different strategies to approximate the
matrix eΔtH.
Taylor expansion: Based on a Taylor expansion of the
matrix eΔtH = ∑k=0,...,∞
(ΔtH)k/k!, we have proposed a novel family of
high-order explicit leap-frog time schemes,
called LFN, where N denotes the order of the
leap-frog scheme. These time schemes ensure the conservation of the
electromagnetic energy as well as the stability under some CFL-like
condition [J4]. The convergence of the semi-discrete
approximation to Maxwell's equations and bounds on the global
divergence error have been rigorously established in
[J3]. This family of time schemes has been
originally proposed for the Maxwell equations in the case of
non-conducting material. Its extension to handle conductive materials
has been studied in [IJCM:10] including the stability analysis.
Chebyshev expansion: A well-konwn alternative for Taylor
expansion is to use Chebyshev polynomials to construct approximations
to Φ(Δt). The basic idea is to expand the time
evolution matrix for a specific time-step Δt in matrix
valued Chebyshev polynomials on the domain of eigenvalues of H,
which lies entirely on the imaginary axis since H is
skew-symmetric. For proper application of the expansion, the domain
of eigenvalues is rescaled to [-1,1], by considering the
matrix B=iH/∥H∥1,
where ∥H∥1 denotes the 1-norm of the
matrix. Operating on state Y(0), the expansion becomes
where I is the identity
matrix, z=Δt∥H∥, Jk is
the kth-order Bessel function and Tk(B) =
ikCk(B). In practice, the summation in
Eq. (2) should be truncated at some expansion
index K. This number depends on the value of z, since
the amplitude of the coefficients Jk(z) decrease
exponentially for k>z. This one-step integrator can be viewed
as an extremely stable time-integration algorithm because it yields an
approximation to the exact time evolution operator that is exact to
nearly machine precision.
Explicit time integration schemes are subjected to stability
conditions that become very restrictive when the underlying mesh is
locally refined since the global time step is deduced from the volume
of the smallest mesh element. Two main strategies can be considered to
improve this situation: local time stepping and implicit time
integration. Although the adoption of an implicit time integration
scheme will allow to overcome the restrictive constraint on the time
step for locally refined meshes, it is still not clear whether the
resulting numerical methodology will demonstrate a superiority in
terms of accuracy and overall computing cost over the original
methodology based on an explicit time integration scheme. On one hand,
the dispersion error should be minimized while taking care to the
increase in complexity of the time integration technique. On the other
hand, the computing cost is also directly impacted by the fact that at
each time step, an implicit time integration scheme yields the
inversion of a large sparse linear system. For certain linear systems
of partial differential equations, the matrix of this system is
constant as far as the time step is fixed during the simulation. The
linear system solver can certainly exploit this fact but this will
probably not translate into a drastic reduction of the cost of a
single time step. Taking into account all these issues, a locally
implicit time integration scheme could be the best compromise when
solving an unsteady wave propagation problem on a locally refined
mesh. The resulting hybrid explicit/implicit time integration strategy
raises several challenges both from the mathematical analysis
viewpoint (stability and accuracy, especially for what concern
numerical dispersion) and from the computer implementation viewpoint
(data structures, parallel computing aspects). Our activities in this
domain aim at the design of hybrid explicit/implicit integration
schemes in conjunction with high order discontinuous Galerkin methods
on locally refined triangular or tetrahedral meshes [J5] - [R5] - [P1] - [C5].
Most discontinuous Galerkin methods developed so far do not include
the study of the error that is due to the discretization of the
geometry (the so-called geometrical error). Actually, the previous
works in DG method consists in mapping, under a linear bijective
transformation, all elements in the physical domain onto a single
reference element for which the local DG matrices are precomputed and
stored once and for all. This technique is subject to certain
constrains on the geometry of physical elements, e.g., straight-sides
for triangles or planar faces for tetrahedra. I have recently
developed a high-order geometrical mapping combined with high-order
DGTD method for the solution of
the Maxwell equations on curvilinear domains. For elements not
intersecting the curved boundary or any curved surface inside the
domain, a standard interpolation and numerical integration are
used. But curved elements are treated through a high-order geometrical
mapping which is based on three ingredients: (a) a high-order
approximation of the curved boundary, (b) a geometric adaptation of
the nodal points inside curved elements, and (c) a proper numerical
integration scheme to evaluate the local DG matrices. I have showed in
[J6] that the DG method based on linear
transformation is inaccurate for curved domains, and
that a higher-order boundary representation introduces a dramatic
improvement in the accuracy of the numerical approximation.
There exist several propagation settings for which the use of a single
geometrical element type (the DGTD method that we have developed so
far is based on a simplex) in the computational domain discretization
process may not be optimal. Instead, one would ideally allow the
combined use of different element types e.g. quadrangles and triangles
in the 2D case, or hexahedra and tetrahedra in the 3D case, possibly
in a non-conforming way (i.e., allowing hanging nodes). Our main
objective is to develop a high-order DGTD-PpQq
method based on non-conforming hybrid hexahedral/tetrahedral meshes
for the numerical simulation of 3D time domain electromagnetic wave
propagation problems.