Transparent boundary conditions for a Discontinuous Galerkin Trefftz method

Transparent boundary conditions in a Discontinuous Galerkin Trefftz method

Abstract

The modeling and simulation of electromagnetic wave propagation is often accompanied by a restriction to bounded domains which requires the introduction of artificial boundaries. The corresponding boundary conditions should be chosen in order to minimize parasitic reflections. In this paper, we investigate a new type of transparent boundary condition for a discontinuous Galerkin Trefftz finite element method. The choice of a particular basis consisting of polynomial plane waves allows us to split the electromagnetic field into components with a well specified direction of propagation. The reflections at the artificial boundaries are then reduced by penalizing components of the field incoming into the space-time domain of interest. We formally introduce this concept, discuss its realization within the discontinuous Galerkin framework, and demonstrate the performance of the resulting approximations by numerical tests. A comparison with first order absorbing boundary conditions, that are frequently used in practice, is made. For a proper choice of basis functions, we observe spectral convergence in our numerical test and an overall dissipative behavior for which we also give some theoretical explanation.

keywords:
transparent boundary conditions, discontinuous Galerkin method, finite element method, Trefftz methods, electrodynamics, wave propagation
1

1 Introduction

We consider the propagation of electromagnetic waves in a domain filled with a non-conducting dielectric medium. In the absence of charges and source currents, the evolution of the electromagnetic fields is governed by the time-dependent Maxwell equations

(1)

The electric permittivity and the magnetic permeability are assumed to be piecewise constant. At time the electric and magnetic fields and are prescribed by the initial conditions

(2)

If the fields satisfy the constraint conditions and in the beginning, then

(3)

which follows by taking the divergence in (1). The two constraint conditions in (3) express the absence of electric charges and magnetic monopoles, respectively. If the computational domain is bounded, the system has to be complemented by appropriate boundary conditions. We will consider different types of conditions that all can be cast in the general abstract form

(4)

here is the outward directed unit normal vector at the domain boundary.

Problems that are described by such a system of equations arise in various applications, for instance, in the modeling of optical wave guides Hadley92 () or in antenna design Milligan05 (). In such cases, boundary conditions of the form

(5)

may be used to model various physically relevant situations, e.g., the presence of perfect electric and magnetic conductors or the action of surface currents describing the emission of energy by an antenna, but also the presence of artificial boundaries resulting from a truncation of the domain which is often introduced to make a simulation feasible. Following the physical intuition, appropriate boundary conditions at such artificial boundaries should allow waves to leave the domain without significant reflection. The first order absorbing boundary condition

(6)

is widely used for this purpose; here is the intrinsic impedance of the medium. This condition mimics the Silver-Müller radiation condition Muller52 (); Barucq97 (); Li14 (), and it is satisfied exactly by plane waves propagating in the outward normal direction. A brief inspection of the Poynting vector

reveals that energy is dissipated by transmission through the boundary at every point on the boundary. We will refer to this condition as first-order absorbing or Silver-Müller condition throughout the paper.

The simple choice (6) can be improved in several ways: In Joly89 (), a more accurate absorbing boundary condition is formulated that still involves only first order derivatives of the fields; for a stability analysis, see also Sonnendrucker10 (). Other possibilities include the classical Bayliss-Turkel and Enquist-Majda conditions Turkel80 (); Turkel85 () and  Engquist77 (); Engquist79 (), which allow to systematically construct conditions for arbitrary order. Due to lack of stability, these are however hardly ever used in practice. Let us also mention more recent approaches developed by Warburton, Hagstrom, Higdon, and others Warburton04 (); Hagstrom07 (); Hagstrom98 (); Higdon86 (); Higdon87 (), the pole condition for the Dirichlet-to-Neumann operator, or the use of infinite elements. Another strategy to minimize reflections from the artificial boundaries is to add an exterior absorbing layer, in which the fields decay very fast. This approach, known as perfectly matched layers, has been used very successfully in practice Berenger94 (); Berenger96 (). The appropriate choice of geometric and physical parameters of the absorbing layer is however not always completely clear in practice, and in some cases it may be necessary to extend the computational domain substantially. In principle, it is also possible to formulate exact boundary conditions, e.g., by the coupling to a boundary integral formulation for the exterior domain Claes80 (); Song87 (). This treatment leads to boundary conditions that are non-local in space and/or time Hiptmair03 (); Kurz99 (), which complicates numerical realization. For a review and a comparison of various kinds of non-absorbing, transparent, or non-reflecting boundary conditions, let us refer to zschiedrichphd (); deaphd () and the references given there.

In this paper, we follow a different strategy for devising local transparent boundary conditions. The intuition behind our approach is the following: Motivated by some of the approaches mentioned above, we assume that at any point of the boundary the electromagnetic fields can be expanded into, or at least approximated by, a superposition of plane waves propagating into specific directions . The three vectors ,, and are assumed to be normalized and orthogonal. If the wave is not reflected at the boundary, one would expect that

(7)

Incorporating such a condition in an adequate manner into a numerical scheme should therefore help to suppress non-physical reflections at artificial boundaries. A similar idea has been used previously in the context of finite difference Trefftz schemes Tsukerman06 (); Tsukerman-book07 (). To evaluate the stability of such a boundary condition, let us again consider the energy flux

across the boundary. Note that because of (7), the summation only runs over indices with and . If the wave at the boundary is mainly propagating in one out-going direction, one can argue that the last term is dominated by the first term on the right hand side, and one obtains outflow of energy over the boundary.

In order to incorporate a boundary condition related to (7) into a numerical method, one has to be able to split the approximation of the electromagnetic field locally into plane waves. This could be realized within the framework of generalized finite elements Babuska92 (); Melenk96 () or via Trefftz finite difference approximations Tsukerman05 (); Tsukerman-book07 (); Tsukerman06 (). Another possibility is provided by the discontinuous Galerkin framework Reed1973 (); Cockburn2001 (); Fezoui2005 (), which allows one to systematically couple almost arbitrary local approximations for the simulation on the global level and to incorporate rather general boundary and interface conditions by some sort of penalization.

In this paper, we consider a space-time discontinuous Galerkin framework for Maxwell’s equations similar to that introduced in Monk14 (); Lilienthal14 (), and we utilize polynomial Trefftz functions for the local approximation which satisfy (1)-(3) exactly on every element. This results in a discontinuous Galerkin Trefftz method that has previously been described in (1+1) dimensions Kretzschmar14 () and later in (3+1) dimensions sisc_2014 (); see also Petersen:2009id (); Farhat14 () for a related Trefftz method in acoustics. The numerical approximation of partial differential equations by Trefftz functions has been proposed in trefftz1926 () and since then been investigated intensively; see e.g. Ruge89 (); Jirousek97 (); Herrera2000 (). Since for the problem under investigation, the Trefftz functions depend on space and time, we automatically arrive at a space-time method. Let us refer to Ruge89 (); Herrera2000 () for a review on the topic and also to Huttunen07 (); Badics08 (); Moiola:2011io () for wave propagation problems in the frequency domain.

One of the basic building blocks of our method is the explicit construction of a basis for the local Trefftz spaces consisting of polynomial plane waves. This allows us to obtain the required local splitting of the discrete electromagnetic fields into plane waves. The second step consists in formulating a variational form of the absorbing boundary condition (7) that can be incorporated within the discontinuous Galerkin framework. Similar to the realization of other boundary conditions in a discontinuous Galerkin method, the condition (7) will be satisfied approximately by some sort of penalization.

To illustrate the benefits of our approach, we present numerical tests including a comparison with first order absorbing boundary conditions. In our computations, we observe spectral convergence for a model problem, provided that the propagation directions of the polynomial plane wave basis functions are chosen appropriately. This indicates that the boundary condition may formally be accurate of arbitrary order. With our numerical tests, we also illustrate energy dissipation and thus stability of the absorbing boundary conditions.

The outline of the paper is as follows: In Section 2, we introduce the space-time discontinuous Galerkin framework which is the basis for our numerical method. In Section 3, we then construct the plane wave basis for the local Trefftz approximation spaces and we sketch the construction for two-dimensional problems underlying our numerical tests. The implementation of the new transparent boundary condition is discussed in detail in Section 4, and results of numerical tests for two simple test problems are reported in Section 5. The presentation closes with a short summary.

2 A space-time Discontinuous Galerkin formulation

For the numerical simulation of the initial boundary value problem (1)–(4), we consider a space-time discontinuous Galerkin method. We utilize Trefftz polynomials for the local approximations which, by definition, satisfy Maxwell’s equations exactly. An appropriate choice of the basis allows us to expand the numerical solution locally into polynomial plane waves and to apply our new transparent boundary condition. In this section, we introduce the general framework of the method. The construction of a plane wave basis for the polynomial Trefftz space and incorporation of the boundary conditions will be addressed in the following two sections.

2.1 Notation

Let be a non-overlapping partition of the domain into regular elements , e.g., tetrahedral, parallelepipeds, prisms, etc. We denote by the set of element interfaces and by the set of faces on the boundary. On an element interface , any piecewise smooth function takes on two values and . We then denote by

the average and the jump of the tangential component of across , respectively; here denotes the outward normal vector on the boundary of the element . Now let be a partition of the time axis into intervals . For every space-time element with , we denote by the space of polynomials in four variables with order up to . We assume that and are constant on , and call

(8)

the space of local Trefftz polynomials; this is the space of vector valued polynomials up to order satisfying Maxwell’s equations (1) and the constraint conditions (3) exactly on the corresponding space-time element.

2.2 The space-time DG framework

For the discretization of the wave propagation problem (1)–(4), we consider a space-time discontinuous Galerkin framework in the spirit of Monk14 (); Lilienthal14 (), but with different approximation spaces and a particular choice of numerical fluxes. On every time slab , we approximate the field by piecewise polynomial Trefftz functions in

(9)

Using these approximation spaces in a space-time discontinuous Galerkin framework of sisc_2014 () yields

Method 1 (Space-time discontinuous Galerkin Trefftz method).


Set , . For find such that for all

This scheme amounts to the methods presented in Kretzschmar14 (); sisc_2014 () with a particular choice of numerical fluxes. Note that, in order to complete the definition of Method 1, we still have to specify the bilinear and linear terms and that account for the boundary conditions. This will be done in Section 4.

2.3 Basic properties of the method

Before we proceed, let us make some general remarks about this numerical scheme; see sisc_2014 () for details and proofs.

(i) Since the approximating functions satisfy Maxwell’s equations exactly on every element, the formulation only contains spatial and temporal interface terms, which penalize the tangential discontinuity of the fields.

(ii) Assume that the true solution of problem (1)–(4) is sufficiently smooth and that the boundary terms are consistently chosen, e.g., such that holds for every point on the boundary. Under this assumption, the whole method is consistent, i.e., any smooth solution of the problem (1)–(4) also satisfies the discrete variational principle. To see this, let us have a closer look onto the discrete variational problem: by tangential continuity of the fields, the last term in the first line drops out. Due to continuity in time, the first terms of the first and third line cancel, whereas the boundary terms cancel by assumption.

(iii) Under mild conditions on the boundary terms, the discrete variational problem for one time slab can be shown to be well-posed. Let denote the spatial mesh-size, the size of the time step, and the spatial dimension. Then the first two terms, which are symmetric positive definite, scale like while the interface and boundary terms scale like . Therefore, the left hand side of the variational principle defines an elliptic bilinear form provided that the time step size is not too large. The smallness condition on can be dropped, if the boundary terms are dissipative in nature, which is the case for many relevant conditions; see the remark at the end of (iv).

(iv) The following energy identity holds

This can be seen from adding a zero term to the variational principle, testing with and , and some elementary algebraic manipulations; see sisc_2014 () for details and proofs. The boundary term with arises from partial integration of one curl operator.

(v) Assume that , which we call condition (D). Then the discrete variational problem is well-posed without any restriction on the time step size. Condition (D) is in fact valid for various types of boundary conditions; see Section 4 for details. If additionally , then the discrete electromagnetic energy defined by is monotonically decreasing in time. We therefore call boundary conditions having the property (D) of dissipative nature.

For details and proofs and some further properties of the resulting scheme, let us refer to  sisc_2014 (); similar results for related discontinuous Galerkin methods based on more standard polynomial spaces can be found in Monk14 (); Lilienthal14 ().

2.4 Implementation

Method 1 yields an implicit time stepping scheme. To evolve the discrete solution from time step to time step , one has to solve a linear system corresponding to the discrete variational problem. Let us sketch the basic structure of this system: After choosing a basis for the piecewise Trefftz space , we can expand the approximate solution with respect to this basis into . The discrete variational problem of Method 1 is then equivalent to the linear system

with matrices , , and vector defined by

According to point (iii) in the discussion of Section 2.3, the matrix is positive definite, provided that the time step is not too large in comparison with the mesh size. For dissipative boundary conditions, this holds without restriction on the size of the time step. In general, will however not be symmetric.

Also note that up to translation of time, the same basis can be used on every time slab. Therefore, if the size of the time step is kept constant, i.e., for all with some , then the matrices and are independent of . This situation is particularly convenient from a computational point of view, since a factorization of the matrix may then be computed once a-priori, and the update of the coefficient vectors from step to only requires one matrix-vector multiplication and a forward-backward substitution. Even if no factorization of is available, the linear system can be solved with acceptable computational effort by some iterative method, as the matrix stems from discretization of a linear hyperbolic problem and therefore usually has a moderate condition number.

3 A basis for the space of Trefftz polynomials

For the local approximation of the electromagnetic fields on every space-time element , we use vector valued polynomials satisfying Maxwell’s equations (1) and the constraint conditions (3) exactly. In this section, we construct a particular basis for this space of Trefftz polynomials consisting of polynomial plane waves, and we discuss some basic properties of this construction. Since we only consider single elements , we assume throughout this section that the material parameters and are positive constants.

3.1 Polynomial plane wave functions

As a basis for the local Trefftz space on the element , we consider polynomial plane waves of the form

(10)

The hat symbols are used to denote constant vectors of unit length. Note that the material properties enter explicitly via the intrinsic impedance and the speed of light . Therefore, the Trefftz basis naturally adapts to local changes in the material.

Lemma 2.

Assume that either or that and , , are mutually orthogonal with . Then any function of the form is a vector valued polynomial of degree which satisfies Maxwell’s equations (1) and the constraint conditions (3). We call a polynomial plane wave.

Proof.

The case yields constant functions and the assertion is clear. Now assume : By definition, the electric field component has the form . One can then verify by direct computation that and ; similar expressions are obtained for the magnetic field component. Maxwell’s equations (1) then reduce to the algebraic conditions

and

The two equations are satisfied if and that which are the assumptions of the Lemma. Additionally, we have . Therefore, also the constraint conditions are satisfied if the directions , , are orthogonal. ∎

3.2 The polynomial Trefftz space

We will now utilize the polynomial plane wave functions introduced in the previous section to define a special basis for the space of local Trefftz polynomials.

Lemma 3.

Let be constant on . Then

  1. .

  2. For there exist orthogonal vector triples , of unit length with and such that the functions , are linearly independent.

  3. The functions , , , form a basis of .

Proof of assertions 1. and 3..

The first assertion follows from the explicit construction of a basis for the space of divergence free Trefftz polynomials in sisc_2014 (). Note that for polynomials of four variables we have . The two Maxwell equations give independent conditions. Applying the divergence operator to Maxwell’s equations yields and . The two constraint conditions thus only have to be required at one point in time and therefore give additional independent conditions. The fact that the two sets of conditions are independent can be seen from the construction in sisc_2014 (), and the assertion follows by counting arguments. Trefftz polynomials with different orders are linearly independent. The third assertion then follows from 1. and 2. by counting the dimensions. ∎

Remark 4.

We cannot provide a complete proof for assertion 2. of Lemma 3 yet. The fact that there exist such linear independent functions can however be verified numerically for any required in practice. An explicit construction, for which we verified the assertion for , is given in Section 3.4 below. Note that by construction and Lemma 2, we know that the polynomials have order and are members of . Moreover, by assertion 1. of the previous lemma, there cannot exist more than linear independent Trefftz polynomials of order .

Remark 5.

According to Lemma 3, every local Trefftz polynomial can be split into a superposition of polynomial plane wave functions . This is the basic requirement for the implementation of the boundary condition (7). If we use the polynomial plane wave basis in the implementation of the method, then the required decomposition is readily available. Let us emphasize that the Trefftz polynomials have coupled electric and magnetic field components, they are functions of space and time, but do not have a tensor-product structure.

3.3 The two-dimensional setting

In cases of translational invariance along one direction, Maxwell’s equations can be cast in a quasi two-dimensional form. For illustration and later reference, let us consider one such case in more detail: We assume that the domain and all fields are homogeneous in the -direction and that the electric field is polarized in this direction. The electromagnetic fields then have the form and with , , and only depending on and . This is known as the TM mode. The computational domain is with being the relevant slice of the three dimensional domain at any fixed and some interval. According to the symmetry assumption on the fields, we require at . Note that this setup still describes a truly three-dimensional problem with symmetry in the direction. We denote by a mesh of the two-dimensional domain , set , and define

(11)

The prime is used here to distinguish this formulation from a fully three-dimensional problem. The construction of a basis for the polynomial Trefftz space is similar to the general case. Here we consider functions of the form

(12)

where and , are orthogonal unit vectors in the - plane. Note that under these assumptions the vector is already fixed by the choice of and the condition .

Lemma 6.

Let be constant on . Then

  1. .

  2. For every there exist orthogonal vector triples , consisting of mutually orthogonal vectors with and . such that the system of functions , are linearly independent.

  3. The functions , , form a basis of .

The proof follows by counting arguments as in the three dimensional case. A particular choice of directions for the two-dimensional setting will again be given in Section 3.4.

Remark 7.

It suffices to consider the field components , , and as functions of , , and only. One could therefore also utilize the alternative representation

(13)

The symbol here denotes the vector-to-scalar and scalar-to-vector curl, respectively. Note that the constraint condition for is satisfied automatically since only depends on and , and the corresponding field points into -direction. The space is isomorphic with and the results stated in the previous lemma carry over.

3.4 Choice of directions

To complete the description of the construction of our basis, we have to find a proper set of independent directions . Let us discuss now in some detail a particular choice that we used to define the polynomial plane wave basis in our numerical experiments.

Three dimensional setting

For we choose six independent constant functions, one for each vector component. For order , we proceed as follows:

  1. We choose distinct numbers , , well distributed in the interval , such that

    This ordering results in an adequate distribution of the directions, cf. Moiola:2011io () and see Figure 1 below.

  2. For every , we choose equidistantly spaced points on the circle . This yields in total different directions , .

  3. For every direction we set and choose two independent, e.g., mutually orthogonal, polarizations , orthogonal to .

  4. For any pair , we finally define .

A possible choice of the directions , is depicted in Figure 1. Note that the linear independence of the corresponding functions can always be verified numerically. A similar construction of directions has been used in moiolaphd () to generate a basis for the time harmonic problem. The analysis of moiolaphd () shows that even (almost) any random choice of directions will yield a linearly independent system of Trefftz functions.

Figure 1: Choice of directions for in the three-dimensional case. The planes corresponding to , are highlighted in different colors. Due to the good distribution and ordering of the levels , the directions are distributed well over the whole sphere.

Two dimensional setting

For the setting described in Section 3.3, the construction of a suitable set of directions is much simpler. We choose three constant functions for , and proceed for as follows:

  1. We choose equidistantly spaced directions , on the unit circle .

  2. We define and .

Linear of the corresponding plane wave functions can again easily be verified numerically.

4 Incorporation of the boundary conditions

To complete the definition of the discontinuous Galerkin method, we now demonstrate how to incorporate various types of boundary conditions. We start by discussing two different implementations for the impedance boundary condition (5), which allow us to treat the perfect-electric-conducting (PEC) and perfect-magnetic-conducting (PMC), as well as the first order absorbing Silver-Müller (SM) boundary condition (6). Our implementation of the transparent boundary condition (7) will turn out to have a very similar structure. In addition to the formulation of these conditions, we also comment on their stability.

4.1 Representation of PEC-like boundary conditions

Let us first consider the impedance boundary condition of the form

(14)

For and , we arrive at the condition for a perfect electric conductor, which is why we call conditions of this form PEC-like. The choice and corresponds to the first-order absorbing boundary condition (6). To incorporate conditions of the form (14) in Method 1, we choose

This form is consistent with the boundary condition (14), i.e., holds if (14) is valid. Let us now consider the energy balance in (iv) of Section 2.3: Testing with yields

For , the last term yields a negative contribution to the energy identity (iv), and the resulting method is stable without restriction on the size of the time step. For and , we thus obtain energy decay on the discrete level.

4.2 Representation of PMC-like boundary conditions

Taking the cross product with from the right in equation (14), it is possible to obtain an alternative equivalent form of the impedance boundary condition, namely

(15)

For and , we arrive at the perfect-magnetic-conducting condition. The choice and yields an equivalent form of the first-order absorbing boundary condition (6). The condition (15) can be incorporated consistently in the discontinuous Galerkin Trefftz method by choosing

Testing with , the boundary term in the energy identity (iv) of Section 2.3 now gives

For any and , we get a negative contribution in the energy identity and thus a dissipative boundary condition. The resulting discrete variational system is then well-posed without restriction on the size of time step.

4.3 A first order absorbing boundary condition

A linear combination of the two conditions (14) and (15) with , , and yields

and , which serves as implementation of the Silver-Müller boundary conditions (6) in our tests.

4.4 New transparent boundary conditions

The basic idea behind our proposal for a transparent boundary condition is to locally expand the electromagnetic field into a superposition of plane waves

(16)

and then suppress the incoming parts by appropriate penalization. Since we are using a basis consisting of plane waves , such a decomposition of the discretized fields is readily available. For the approximation of the transparent boundary condition (7) within our discontinuous Galerkin framework, we then consider the choice

(17)

and we set . Here denotes the incoming part of the electromagnetic fields, i.e., summation is done only over indices with . This condition has a similar form as the Silver-Müller condition stated above, but only the incoming fields are taken into account. Let us also examine the energy balance for the new boundary condition: Testing with , and assuming for simplicity, we obtain

Here, and denote the out-going field components. Similar as on the continuous level, we may argue that the first two terms in the second line will give a positive contribution, if the numerical solution is mainly directed into an outward direction. This can be expected to be the case, if the continuous solution has this behavior. The third and fourth term can then be absorbed into the first two and the last term via a Young’s inequality. In summary, we thus expect a decay of the discrete energy, which is what we actually observe in our numerical tests.

5 Test problems and numerical results

For numerical validation of the new transparent boundary condition, we consider two test problems. The first test problem studies the propagation of a plane wave. In this scenario an analytic solution is available, which allows us to conduct a numerical convergence study. In the second test problem, we consider the propagation of a cylindrical wave. We evaluate the effect of the transparent boundary condition on the dissipation of energy by comparing the numerical solutions obtained on a large domain and on an artificially truncated domain with different choices of boundary conditions.

5.1 Transmission of a Plane Wave

We consider a plane wave propagating in direction through a homogeneous medium with parameters . The fields and with

(18)

satisfy Maxwell’s equations (1) and the constraint conditions (3), and they will serve as the reference solution. The evolution of the field component of the analytic solution over time is depicted in Fig. 2.

Figure 2: Evolution of the electric field component of the plane wave (18) propagating through the computational domain .

Since the fields and are independent of the third coordinate direction, it suffices to consider a geometrically two-dimensional setting; see Section 3.3. As a computational domain, we choose . On the incoming boundaries at and , the fields are set to that of the analytic solution by the PEC-like boundary condition (14) with and determined from the analytic solution. Different kinds of boundary conditions are utilized at the boundaries and , where the wave leaves the domain.

For our simulations, we start from a uniform initial mesh with a mesh size resulting in rectangular elements. The size of the time step is chosen as throughout our tests. We employ Method 1 with approximation spaces and different choices of . According to our considerations in Section 3.3, the total number of degrees of freedom for one time step is then . Simulations are carried out until , where the wave should have left the domain almost completely.

In a first series of tests, we evaluate the order of convergence with respect to refinement of the spatial and temporal mesh size. We run simulations for different approximation orders and on sequences of uniformly refined meshes. In all tests, is utilized as the time step. In Fig. 3 we display the relative error of the computed approximations in the space-time norm as a function of the mesh size .

Figure 3: Relative error as a function of the spatial mesh-size for orders .

In all simulations, we observe convergence rates of order which are optimal with respect to the approximation properties of a piecewise polynomial space of order . The discontinuous Galerkin Trefftz method thus yields a quasi optimal approximation.

To evaluate more closely the effect of the transparent boundary condition, we display in Fig. 4 the errors obtained for different choices of boundary conditions. In particular, we compare the implementation of the Silver-Müller condition given at the end of Section 4.2 with the transparent boundary condition discussed in Section 4.4. Simulations are carried out with polynomial approximation orders , and on a uniform mesh with mesh-size , , and elements.

Figure 4: Relative error