Stabilization in relation to wavenumber in HDG methods
SSTABILIZATION IN RELATION TO WAVENUMBER IN HDGMETHODS
J. GOPALAKRISHNAN, S. LANTERI, N. OLIVARES, AND R. PERRUSSEL
Abstract.
Simulation of wave propagation through complex media relies on properunderstanding of the properties of numerical methods when the wavenumber is realand complex. Numerical methods of the Hybrid Discontinuous Galerkin (HDG) typeare considered for simulating waves that satisfy the Helmholtz and Maxwell equations.It is shown that these methods, when wrongly used, give rise to singular systems forcomplex wavenumbers. A sufficient condition on the HDG stabilization parameter forguaranteeing unique solvability of the numerical HDG system, both for Helmholtz andMaxwell systems, is obtained for complex wavenumbers. For real wavenumbers, resultsfrom a dispersion analysis are presented. An asymptotic expansion of the dispersionrelation, as the number of mesh elements per wave increase, reveal that some choices ofthe stabilization parameter are better than others. To summarize the findings, thereare values of the HDG stabilization parameter that will cause the HDG method tofail for complex wavenumbers. However, this failure is remedied if the real part of thestabilization parameter has the opposite sign of the imaginary part of the wavenumber.When the wavenumber is real, values of the stabilization parameter that asymptoti-cally minimize the HDG wavenumber errors are found on the imaginary axis. Finally,a dispersion analysis of the mixed hybrid Raviart-Thomas method showed that itswavenumber errors are an order smaller than those of the HDG method.
Background
Wave propagation through complex structures, composed of both propagating andabsorbing media, are routinely simulated using numerical methods. Among the vari-ous numerical methods used, the Hybrid Discontinuous Galerkin (HDG) method hasemerged as an attractive choice for such simulations. The easy passage to high orderusing interface unknowns, condensation of all interior variables, availability of error es-timators and adaptive algorithms, are some of the reasons for the adoption of HDGmethods.It is important to design numerical methods that remain stable as the wavenumbervaries in the complex plane. For example, in applications like computational lithography,one finds absorbing materials with complex refractive index in parts of the domain of
Key words and phrases.
HDG, Raviart-Thomas, dispersion, dissipation, absorbing material, com-plex, wave speed, optimal, stabilization, Helmholtz, Maxwell, unisolvency.JG and NO were supported in part by the NSF grant DMS-1318916 and the AFOSR grant FA9550-12-1-0484. NO gratefully acknowledges support in the form of an INRIA internship where discussionsleading to this work originated. All authors wish to thank INRIA Sophia Antipolis M´editerran´ee forhosting the authors there and facilitating this research. a r X i v : . [ m a t h . NA ] M a r J. GOPALAKRISHNAN, S. LANTERI, N. OLIVARES, AND R. PERRUSSEL simulation. Other examples are furnished by meta-materials. A separate and importantreason for requiring such stability emerges in the computation of resonances by iterativesearches in the complex plane. It is common for such iterative algorithms to solve asource problem with a complex wavenumber as its current iterate. Within such algo-rithms, if the HDG method is used for discretizing the source problem, it is imperativethat the method remains stable for all complex wavenumbers.One focus of this study is on complex wavenumber cases in acoustics and electromag-netics, motivated by the above-mentioned examples. Ever since the invention of theHDG method in [4], it has been further developed and extended to other problems inmany works (so many so that it is now impractical to list all references on the subjecthere). Of particular interest to us are works that applied HDG ideas to wave propaga-tion problems such as [6, 8, 10, 11, 12, 13, 14]. We will make detailed comparisons withsome of these works in a later section. However, none of these references address thestability issues for complex wavenumber cases. While the choice of the HDG stabiliza-tion parameter in the real wave number case can be safely modeled after the well-knownchoices for elliptic problems [5], the complex wave number case is essentially different.This will be clear right away from a few elementary calculations in the next section,which show that the standard prescriptions of stabilization parameters are not alwaysappropriate for the complex wave number case. This then raises further questions onhow the HDG stabilization parameter should be chosen in relation to the wavenumber,which are addressed in later sections.Another focus of this study is on the difference in speeds of the computed and theexact wave, in the case of real wavenumbers. By means of a dispersion analysis, one cancompute the discrete wavenumber of a wave-like solution computed by the HDG method,for any given exact wavenumber. An extensive bibliography on dispersion analysesfor the standard finite element method can be obtained from [1, 7]. For nonstandardfinite element methods however, dispersion analysis is not so common [9], and for theHDG method, it does not yet exist. We will show that useful insights into the HDGmethod can be obtained by a dispersion analysis. In multiple dimensions, the discretewavenumber depends on the propagation angle. Analytic computation of the dispersionrelation is feasible in the lowest order case. We are thus able to study the influenceof the stabilization parameter on the discrete wavenumber and offer recommendationson choosing good stabilization parameters. The optimal stabilization parameter valuesare found not to depend on the wavenumber. In the higher order case, since analyticcalculations pose difficulties, we conduct a dispersion analysis numerically.We begin, in the next section, by describing the HDG methods. We set the stage forthis study by showing that the commonly chosen HDG stabilization parameter valuesfor elliptic problems are not appropriate for all complex wavenumbers. In the subse-quent section, we discover a constraint on the stabilization parameter, dependent on thewavenumber, that guarantees unique solvability of both the global and the local HDGproblems. Afterward, we perform a dispersion analysis for both the HDG method anda mixed method and discuss the results.
TABILIZATION AND WAVENUMBER IN HDG 3
Methods of the HDG type
We borrow the basic methodology for constructing HDG methods from [4] and apply itto the time-harmonic Helmholtz and Maxwell equations (written as first order systems).While doing so, we set up the notations used throughout, compare the formulationwe use with other existing works, and show that for complex wavenumbers there arestabilization parameters that will cause the HDG method to fail.
Undesirable HDG stabilization parameters for the Helmholtz system.
Webegin by considering the lowest order HDG system for Helmholtz equation. Let k be acomplex number. Consider the Helmholtz system on Ω ⊂ R with homogeneous Dirichletboundary conditions, ˆ ık ⃗U + ⃗∇ Φ = ⃗ , in Ω , (1a)ˆ ıkΦ + ⃗∇⋅ ⃗U = f, in Ω , (1b) Φ = , on ∂ Ω , (1c)where f ∈ L ( Ω ) . Let T h denote a square or triangular mesh of disjoint elements K , soΩ = ∪ K ∈T h K , and let F h denote the collection of edges. The HDG method produces anapproximation (⃗ u, φ, ˆ φ ) to the exact solution ( ⃗U , Φ, ˆ Φ ) , where ˆ Φ denotes the trace of Φ on the collection of element boundaries ∂ T h . The HDG solution (⃗ u, φ, ˆ φ ) is in the finitedimensional space V h × W h × M h defined by V h = {⃗ v ∈ ( L ( Ω )) ∶ ⃗ v ∣ K ∈ V ( K ) , ∀ K ∈ T h } W h = { ψ ∈ L ( Ω ) ∶ ψ ∣ K ∈ W ( K ) , ∀ K ∈ T h } M h = { ˆ ψ ∈ L ( ⋃ F ∈F h F ) ∶ ˆ ψ ∣ F ∈ M ( F ) , ∀ F ∈ F h and ˆ ψ ∣ ∂ Ω = } , with polynomial spaces V ( K ) , W ( K ) , and M ( F ) specified differently depending onelement type: Triangles Squares V ( K ) = (P p ( K )) V ( K ) = (Q p ( K )) W ( K ) = P p ( K ) W ( K ) = Q p ( K ) M ( F ) = P p ( F ) M ( F ) = P p ( F ) . Here, for a given domain D , P p ( D ) denotes polynomials of degree at most p , and Q p ( D ) denotes polynomials of degree at most p in each variable.The HDG solution satisfies ∑ K ∈T h ˆ ık (⃗ u, ⃗ v ) K − ( φ, ⃗∇⋅ ⃗ v ) K + ⟨ ˆ φ, ⃗ v ⋅ ⃗ n ⟩ ∂K = , (2a) ∑ K ∈T h −( ⃗∇⋅ ⃗ u, ψ ) K + ⟨ τ ˆ φ, ψ ⟩ ∂K − ⟨ τ φ, ψ ⟩ ∂K − ˆ ık ( φ, ψ ) K = −( f, ψ ) Ω , (2b) ∑ K ∈T h ⟨⃗ u ⋅ ⃗ n + τ ( φ − ˆ φ ) , ˆ ψ ⟩ ∂K = , (2c) J. GOPALAKRISHNAN, S. LANTERI, N. OLIVARES, AND R. PERRUSSEL
Reference Their notations and equations Connection to our formulation[6]Helmholtz case ⃗ q [6] + ⃗∇ u [6] = ⃗ ⃗∇⋅ ⃗ q [6] − k u [6] = q [6] ⋅ ⃗ n = ⃗ q [6] ⋅ ⃗ n + ˆ ıτ [6] ( u [6] − ˆ u [6] ) τ [6] = k τ ˆ ıku [6] = φ ⃗ q [6] = ⃗ u [10]Helmholtz case ˆ ık ⃗ q [10] + ⃗∇ u [10] = ⃗ ıku [10] + ⃗∇⋅ ⃗ q [10] = q [10] ⋅ ⃗ n = ⃗ q [10] ⋅ ⃗ n + τ [10] ( u [10] − ˆ u [10] ) τ [10] = τu [10] = φ ⃗ q [10] = ⃗ u [12]2D Maxwell case ˆ ıω [12] ε r E [12] − ∇× ⃗ H [12] = ıω [12] µ r ⃗ H [12] + ⃗∇× E [12] = ⃗ H [12] = ⃗ H [12] + τ [12] ( E [12] − ˆ E [12] )⃗ t τ [12] = √ ε r µ r τω [12] = ω √ ε µ E [12] = √ ε r E, ⃗ H [12] = √ µ r ⃗ H [14]Maxwell case µ ⃗ w [14] − ⃗∇× ⃗ u [14] = ⃗ ⃗∇× ⃗ w [14] − εω ⃗ u [14] = ⃗ w [14] = ⃗ w [14] + τ [14] (⃗ u [14] − ˆ u [14] ) × ⃗ n τ [14] = ˆ ı √ εω µ τµ ⃗ w [14] = − ˆ ık ⃗ H, with k = ω √ µε, ⃗ u [14] = ⃗ E Table 1.
Comparison with some HDG formulations in other papers. No-tations in the indicated external references are used after subscriptingthem by the reference number. Notations without subscripts are thosedefined in this paper.for all ⃗ v ∈ V h , ψ ∈ W h , and ˆ ψ ∈ M h . The last equation enforces the conservativity of thenumerical flux ˆ u ⋅ ⃗ n = ⃗ u ⋅ ⃗ n + τ ( φ − ˆ φ ) . (3)The stabilization parameter τ is assumed to be constant on each ∂K . We are interestedin how the choice of τ in relation to k affects the method, especially when k is complexvalued. Comparisons of this formulation with other HDG formulations for Helmholtzequations in the literature are summarized in Table 1.One of the main reasons to use an HDG method is that all interior unknowns ( ⃗ u, φ )can be eliminated to get a global system for solely the interface unknowns ( ˆ φ ). This is TABILIZATION AND WAVENUMBER IN HDG 5 possible whenever the local systemˆ ık (⃗ u, ⃗ v ) K − ( φ, ⃗∇⋅ ⃗ v ) K = −⟨ ˆ φ, ⃗ v ⋅ ⃗ n ⟩ ∂K , ∀⃗ v ∈ V ( K ) , (4a) −( ⃗∇⋅ ⃗ u, ψ ) K − ⟨ τ φ, ψ ⟩ ∂K − ˆ ık ( φ, ψ ) K = −⟨ τ ˆ φ, ψ ⟩ ∂K − ( f, ψ ) K , ∀ ψ ∈ W ( K ) , (4b)is uniquely solvable. (For details on this elimination and other perspectives on HDGmethods, see [4].) In the lowest order ( p =
0) case, on a square element K of side length h , if we use a basis in the following order ⃗ u = [ ] , ⃗ u = [ ] , φ = , on K, then the element matrix for the system (4) is M = ⎡⎢⎢⎢⎢⎢⎣ ˆ ık h ık h
00 0 − hτ − ˆ ık h ⎤⎥⎥⎥⎥⎥⎦ . This shows that if 4 τ = − ˆ ıkh, (5)then M is singular, and so the HDG method will fail. The usual recipe of choosing τ = is therefore inappropriate when k is complex valued. Intermediate case of the 2D Maxwell system.
It is an interesting exercise toconsider the 2D Maxwell system before going to the full 3D case. In fact, the HDGmethod for the 2D Maxwell system can be determined from the HDG method for the2D Helmholtz system. The 2D Maxwell system isˆ ık E − ∇× ⃗H = − J, (6a)ˆ ık ⃗H + ⃗∇× E = , (6b)where J ∈ L ( Ω ) , and the scalar curl ∇× ⋅ and the vector curl ⃗∇× ⋅ are defined by ∇× ⃗H = ∂ H − ∂ H = ⃗∇⋅ R ( ⃗H) , ⃗∇× E = ( ∂ E , − ∂ E ) = R ( ⃗∇E ) . Here R ( v , v ) = ( v , − v ) is the operator that rotates vectors counterclockwise by + π / ⃗ r = − R ⃗H , then (6) becomesˆ ık E + ⃗∇⋅ ⃗ r = − J, − ˆ ık ⃗ r + RR ⃗∇E = , which, since RR ⃗ v = −⃗ v (rotation by π ), coincides with (1) with Φ = E , ⃗U = ⃗ r , and f = − J .This also shows that the HDG method for Helmholtz equation should yield an HDGmethod for the 2D Maxwell system. We thus conclude that there exist stabilizationparameters that will cause the HDG system for 2D Maxwell system to fail. J. GOPALAKRISHNAN, S. LANTERI, N. OLIVARES, AND R. PERRUSSEL
To examine this 2D HDG method, if we let ⃗ H and E denote the HDG approximationsfor R ⃗ r and E , respectively, then the HDG system (2) with ⃗ u and φ replaced by − R ⃗ H and E , respectively, gives ∑ K ∈T h −( E, ⃗∇× ⃗ w ) K + ⟨ ˆ E, ⃗ n × ⃗ w ⟩ ∂K − ˆ ık ( ⃗ H, ⃗ w ) K = , ∑ K ∈T h ˆ ık ( E, ψ ) K − ( ⃗∇× ⃗ H, ψ ) K + ⟨ τ ( E − ˆ E ) , ψ ⟩ ∂K = −( ⃗ J , ψ ) Ω , ∑ K ∈T h ⟨ ̂ R ⃗ H ⋅ ⃗ n, ˆ ψ ⟩ ∂K = , for all ⃗ w ∈ R ( V h ) , ψ ∈ W h and ˆ ψ ∈ M h . We have used the fact that −( R ⃗ H ) ⋅ ⃗ n = ⃗ H ⋅ ⃗ t ,where ⃗ t = R ⃗ n the tangent vector, and we have used the 2D cross product defined by ⃗ v × ⃗ n = ⃗ v ⋅ ⃗ t . In particular, the numerical flux prescription (3) implies − ̂ R ⃗ H ⋅ ⃗ n = − R ⃗ H ⋅ ⃗ n + τ ( E − ˆ E ) , where ̂ R ⃗ H denotes the numerical trace of R ⃗ H . We rewrite this in terms of ⃗ H and E , toobtain ˆ H ⋅ ⃗ t = ⃗ H ⋅ ⃗ t + τ ( E − ˆ E ) . One may rewrite this again, as ˆ H × ⃗ n = ⃗ H × ⃗ n + τ ( E − ˆ E ) . (7)This expression is notable because it will help us consistently transition the numericalflux prescription from the Helmholtz to the full 3D Maxwell case discussed next. Acomparison of this formula with those in the existing literature is included in Table 1. The 3D Maxwell system.
Consider the 3D Maxwell system on Ω ⊂ R with a perfectelectrically conducting boundary condition,ˆ ık ⃗E − ⃗∇× ⃗H = − ⃗ J , in Ω , (8a)ˆ ık ⃗H + ⃗∇× ⃗E = ⃗ , in Ω , (8b) ⃗ n × ⃗E = ⃗ , on ∂ Ω , (8c)where ⃗ J ∈ ( L ( Ω )) . For this problem, T h denotes a cubic or tetrahedral mesh, and F h denotes the collection of mesh faces. The HDG method approximates the exact solution ( ⃗E , ⃗H , ˆ E ) by the discrete solution ( ⃗ E, ⃗ H, ˆ E ) ∈ Y h × Y h × N h . The discrete spaces are definedby Y h = {⃗ v ∈ ( L ( Ω )) ∶ ⃗ v ∣ K ∈ Y ( K ) , ∀ K ∈ T h } N h = { ˆ η ∈ ( L (F h )) ∶ ˆ η ∣ F ∈ N ( F ) , ∀ F ∈ F h and ˆ η ∣ ∂ Ω = ⃗ } , TABILIZATION AND WAVENUMBER IN HDG 7 with polynomial spaces Y ( K ) and N ( F ) specified by:Tetrahedra Cubes Y ( K ) = (P p ( K )) Y ( K ) = (Q p ( K )) N ( F ) = { ˆ η ∈ (P p ( F )) ∶ ˆ η ⋅ ⃗ n = } N ( F ) = { ˆ η ∈ (Q p ( F )) ∶ ˆ η ⋅ ⃗ n = } Our HDG method for (8) is ∑ K ∈T h ˆ ık ( ⃗ E, ⃗ v ) K − ( ⃗∇× ⃗ H, ⃗ v ) K + ⟨( ˆ H − H ) × ⃗ n, ⃗ v ⟩ ∂K = −( ⃗ J , ⃗ v ) Ω , ∀⃗ v ∈ Y h , ∑ K ∈T h −( ⃗ E, ⃗∇× ⃗ w ) K + ⟨ ˆ E, ⃗ n × ⃗ w ⟩ ∂K − ˆ ık ( ⃗ H, ⃗ w ) K = , ∀ ⃗ w ∈ Y h , ∑ K ∈T h ⟨ ˆ H × ⃗ n, ˆ w ⟩ ∂K = , ∀ ˆ w ∈ N h , where, in analogy with (7), we now set numerical flux byˆ H × ⃗ n = ⃗ H × ⃗ n + τ ( ⃗ E − ˆ E ) t , (9)where ( ⃗ E − ˆ E ) t denotes the tangential component, or equivalentlyˆ H × ⃗ n = ⃗ H × ⃗ n + τ (⃗ n × ( ⃗ E − ˆ E )) × ⃗ n. Note that the 2D system (6) is obtained from the 3D Maxwell system (8) by assumingsymmetry in x -direction. Hence, for consistency between 2D and 3D formulations, weshould have the same form for the numerical flux prescriptions in 2D and 3D.The HDG method is then equivalently written as ∑ K ∈T h ˆ ık ( ⃗ E, ⃗ v ) K − ( ⃗∇× ⃗ H, ⃗ v ) K + ⟨ τ ( ⃗ E − ˆ E ) × ⃗ n, ⃗ v × ⃗ n ⟩ ∂K = −( ⃗ J , ⃗ v ) Ω , (10a) ∑ K ∈T h −( ⃗ E, ⃗∇× ⃗ w ) K + ⟨ ˆ E, ⃗ n × ⃗ w ⟩ ∂K − ˆ ık ( ⃗ H, ⃗ w ) K = , (10b) ∑ K ∈T h ⟨ ⃗ H + τ ⃗ n × ( ⃗ E − ˆ E ) , ˆ w × ⃗ n ⟩ ∂K = , (10c)for all ⃗ v, ⃗ w ∈ Y h , and ˆ w ∈ N h . For comparison with other existing formulations, seeTable 1.Again, let us look at the solvability of the local element problem ˆ ık ( ⃗ E, ⃗ v ) K − ( ⃗∇× ⃗ H, ⃗ v ) K + ⟨ τ ⃗ E × ⃗ n, ⃗ v × ⃗ n ⟩ ∂K = ⟨ τ ˆ E × ⃗ n, ⃗ v × ⃗ n ⟩ ∂K − ( ⃗ J , ⃗ v ) K , (11a) −( ⃗ E, ⃗∇× ⃗ w ) K − ˆ ık ( ⃗ H, ⃗ w ) K = −⟨ ˆ E, ⃗ n × ⃗ w ⟩ ∂K , (11b)for all ⃗ v, ⃗ w ∈ Y ( K ) . In the lowest order ( p =
0) case, on a cube element K of side length h , if we use a basis in the following order ⃗ E = ⎡⎢⎢⎢⎢⎢⎣ ⎤⎥⎥⎥⎥⎥⎦ , ⃗ E = ⎡⎢⎢⎢⎢⎢⎣ ⎤⎥⎥⎥⎥⎥⎦ , ⃗ E = ⎡⎢⎢⎢⎢⎢⎣ ⎤⎥⎥⎥⎥⎥⎦ , ⃗ H = ⎡⎢⎢⎢⎢⎢⎣ ⎤⎥⎥⎥⎥⎥⎦ , ⃗ H = ⎡⎢⎢⎢⎢⎢⎣ ⎤⎥⎥⎥⎥⎥⎦ , ⃗ H = ⎡⎢⎢⎢⎢⎢⎣ ⎤⎥⎥⎥⎥⎥⎦ , (12) J. GOPALAKRISHNAN, S. LANTERI, N. OLIVARES, AND R. PERRUSSEL then the 6 × M = [( h τ + ˆ ıkh ) I −( ˆ ıkh ) I ] , where I denotes the 3 × τ = − ˆ ıkh, (13)then the local static condensation required in the HDG method will fail in the Maxwellcase also. Behavior on tetrahedral meshes.
For the lowest order ( p =
0) case on a tetrahe-dral element, just as for the cube element described above, there are bad stabilizationparameter values. Consider, for example, the tetrahedral element of size h defined by K = {⃗ x ∈ R ∶ x j ⩾ ∀ j, x + x + x ⩽ h } , (14)with a basis ordered as in (12). The element matrix for the system (11) is then M = ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣( √ + ) h τ + ˆ ıkh −√ h τ −√ h τ −√ h τ ( √ + ) h τ + ˆ ıkh −√ h τ −√ h τ −√ h τ h τ + ˆ ıkh − ˆ ıkh − ˆ ıkh
00 0 0 0 0 − ˆ ıkh ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ . We immediately see that the rows become linearly dependent if ( √ + ) τ = − ˆ ıkh. Hence, for τ = − ˆ ıkh /( √ + ) , the HDG method will fail on tetrahedral meshes.For orders p ⩾
1, the element matrices are too complex to find bad parameter valuesso simply. Instead, we experiment numerically. Setting τ = − ˆ ı , which is equivalent tothe choice made in [14] (see Table 1), we compute the smallest singular value of theelement matrix M (the matrix of the left hand side of (11) with K set by (14)) for arange of normalized wavenumbers kh . Figures 1a and 1b show that, for orders p = p =
2, there are values of kh for which τ = − ˆ ı results in a singular value very close tozero. Taking a closer look at the first nonzero local minimum in Figure 1a, we find thatthe local matrix corresponding to normalized wavenumber kh ≈ .
49 has an estimatedcondition number exceeding 3 . × , i.e., for all practical purposes, the element matrixis singular. To illustrate how a different choice of stabilization parameter τ can affectthe conditioning of the element matrix, Figures 1c and 1d show the smallest singularvalues for the same range of kh , but with τ =
1. Clearly the latter choice of τ is betterthan the former.From another perspective, Figure 1e shows the smallest singular value of the elementmatrix as τ is varied in the complex plane, while fixing kh to 1. Figure 1f is similarexcept that we fixed kh to the value discussed above, approximately 7 .
49. In both cases,
TABILIZATION AND WAVENUMBER IN HDG 9 s m a ll e s t s i ngu l a r v a l ueo f e l e m en t m a t r i x (a) p = τ = − ˆ ı s m a ll e s t s i ngu l a r v a l ueo f e l e m en t m a t r i x (b) p = τ = − ˆ ı s m a ll e s t s i ngu l a r v a l ueo f e l e m en t m a t r i x (c) p = τ = s m a ll e s t s i ngu l a r v a l ueo f e l e m en t m a t r i x (d) p = τ = −20 −10 0 10 20−20i−10i 0i 10i 20i 0.0020.0040.0060.008 (e) kh = p = −4 −2 0 2 4−6i−4i−2i 0i 2i 4i 6i 0.010.020.030.040.050.06 (f) kh ≈ . p = −20 −10 0 10 20−20i−10i 0i 10i 20i 0.0020.0040.0060.0080.01 (g) kh = + ˆ ı , p = −4 −2 0 2 4−6i−4i−2i 0i 2i 4i 6i 0.020.040.060.08 (h) kh ≈ . ( + ˆ ı ) , p = Figure 1.
The smallest singular values of a tetrahedral HDG element matrix we find that the worst values of τ are along the imaginary axis. Finally, in Figures 1gand 1h, we see the effects of multiplying these real values of kh by 1 + ˆ ı . The regionof the complex plane where bad values of τ are found changes significantly when kh iscomplex. Results on unisolvent stabilization
We now turn to the question of how we can choose a value for the stabilizationparameter τ that will guarantee that the local matrices are not singular. The answer,given by a condition on τ , surprisingly also guarantees that the global condensed HDGmatrix is nonsingular. These results are based on a tenuous stability inherited fromthe fact nonzero polynomials are never waves, stated precisely in the ensuing lemma.Then we give the condition on τ that guarantees unisolvency, and before concluding thesection, present some caveats on relying solely on this tenuous stability.As is standard in all HDG methods, the unique solvability of the element problemallows the formulation of a condensed global problem that involves only the interfaceunknowns. We introduce the following notation to describe the condensed systems.First, for Maxwell’s equations, for any η ∈ N h , let ( ⃗ E η , ⃗ H η ) ∈ Y h × Y h denote the fieldssuch that, for each K ∈ T h , the pair ( ⃗ E η ∣ K , ⃗ H η ∣ K ) satisfies the local problem (11) withdata η ∣ ∂K . That is,ˆ ık ( ⃗ E η , ⃗ v ) K − ( ⃗∇× ⃗ H η , ⃗ v ) K + ⟨ τ ⃗ E η × ⃗ n, ⃗ v × ⃗ n ⟩ ∂K = ⟨ τ η × ⃗ n, ⃗ v × ⃗ n ⟩ ∂K , (15a) −( ⃗ E η , ⃗∇× ⃗ w ) K − ˆ ık ( ⃗ H η , ⃗ w ) K = −⟨ η, ⃗ n × ⃗ w ⟩ ∂K , (15b)for all ⃗ v ∈ Y ( K ) , ⃗ w ∈ Y ( K ) . If all the sources in (10) vanish, then the condensed globalproblem for ˆ E ∈ N h takes the form a ( ˆ E, η ) = , ∀ η ∈ N h , (16)where a ( Λ , η ) = ∑ K ∈T h ⟨ ˆ H Λ × ⃗ n, η ⟩ ∂K . By following a standard procedure [4] we can express a (⋅ , ⋅) explicitly as follows: a ( Λ , η ) = ∑ K ∈T h ⟨ ⃗ H Λ × ⃗ n, η ⟩ ∂K + ⟨( ˆ H Λ − ⃗ H Λ ) × ⃗ n, η ⟩ ∂K = ∑ K ∈T h ˆ ık ( ⃗ H Λ , ⃗ H η ) K − ( ⃗∇× ⃗ H Λ , ⃗ E η ) K + ⟨ τ ⃗ n × (⃗ n × ( Λ − ⃗ E Λ )) , η ⟩ ∂K = ∑ K ∈T h ˆ ık ( ⃗ H Λ , ⃗ H η ) K − ˆ ık ( ⃗ E Λ , ⃗ E η ) K + ⟨ τ ⃗ n × ( Λ − ⃗ E Λ ) , ⃗ n × ⃗ E η ⟩ ∂K − ⟨ τ ⃗ n × ( Λ − ⃗ E Λ ) , ⃗ n × η ⟩ ∂K = ∑ K ∈T h ˆ ık ( ⃗ H Λ , ⃗ H η ) K − ˆ ık ( ⃗ E Λ , ⃗ E η ) K − τ ⟨⃗ n × ( Λ − ⃗ E Λ ) , ⃗ n × ( η − ⃗ E η )⟩ ∂K . Here we have used the complex conjugate of (15b) with ⃗ w = ⃗ H Λ , along with the definitionof ˆ H Λ , and then used (15a). TABILIZATION AND WAVENUMBER IN HDG 11
Similarly, for the Helmholtz equation, let (⃗ u η , φ η ) ∈ V h × W h denote the fields suchthat, for all K ∈ T h , the functions (⃗ u η ∣ K , φ η ∣ K ) solve the element problem (4) for givendata ˆ φ = η . If the sources in (2) vanish, then the condensed global problem for ˆ φ ∈ M h iswritten as b ( ˆ φ, η ) = , ∀ η ∈ M h , (17)where the form is found, as before, by the standard procedure: b ( Λ , η ) = ∑ K ∈T h ⟨ ˆ u Λ ⋅ ⃗ n, η ⟩ ∂K = ∑ K ∈T h ˆ ık (⃗ u Λ , ⃗ u η ) K − ˆ ık ( φ Λ , φ η ) K − τ ⟨ Λ − φ Λ , η − φ η ⟩ ∂K . The sesquilinear forms a (⋅ , ⋅) and b (⋅ , ⋅) are used in the main result, which gives suffi-cient conditions for the solvability of the local problems (11), (4) and the global prob-lems (16), (17).Before proceeding to the main result, we give a simple lemma, which roughly speaking,says that nontrivial harmonic waves are not polynomials . Lemma 1.
Let p ⩾ be an integer, ≠ k ∈ C , and D an open set. Then, there is nonontrivial ⃗ E ∈ (P p ( D )) satisfying ⃗∇×( ⃗∇× ⃗ E ) − k ⃗ E = and there is no nontrivial φ ∈ P p ( D ) satisfying ∆ φ + k φ = . Proof.
We use a contradiction argument. If E /≡ ⃗
0, then we may assume without loss ofgenerality that at least one of the components of ⃗ E is a polynomial of degree exactly p .But this contradicts k ⃗ E = ⃗∇×( ⃗∇× ⃗ E ) because all components of ⃗∇×( ⃗∇× ⃗ E ) are polyno-mials of degree at most p −
2. Hence ⃗ E ≡ ⃗
0. An analogous argument can be used for theHelmholtz case as well. (cid:3)
Theorem 1.
Suppose Re ( τ ) ≠ , whenever Im ( k ) = , and (18a)Im ( k ) Re ( τ ) ⩽ , whenever Im ( k ) ≠ . (18b) Then, in the Maxwell case, the local element problem (11) and the condensed globalproblem (16) are both unisolvent. Under the same condition, in the Helmholtz case, thelocal element problem (4) and the condensed global problem (17) are also unisolvent.Proof.
We first prove the theorem for the local problem for Maxwell’s equations. As-sume (18) holds and set ˆ E = ⃗ ⃗ E and ⃗ H must equal ⃗
0. Choosing ⃗ v = ⃗ E , and ⃗ w = ⃗ H , then subtract-ing (11b) from (11a), we getˆ ık (∣∣ ⃗ E ∣∣ K + ∣∣ ⃗ H ∣∣ K ) + ı Im ( ⃗ E, ⃗∇× ⃗ H ) K + τ ∣∣ ⃗ E × ⃗ n ∣∣ ∂K = , whose real part is − Im ( k ) (∣∣ ⃗ E ∣∣ K + ∣∣ ⃗ H ∣∣ K ) + Re ( τ )∣∣ ⃗ E × ⃗ n ∣∣ ∂K = . Under condition (18b), we immediately have that the fields ⃗ E and ⃗ H are zero on K .Otherwise, (18a) implies ⃗ E × ⃗ n ∣ ∂K =
0, and then (11) givesˆ ık ⃗ E − ⃗∇× ⃗ H = , ˆ ık ⃗ H + ⃗∇× ⃗ E = , implying ⃗∇×( ⃗∇× ⃗ E ) = k ⃗ E. By Lemma 1 this equation has no nontrivial solutions in the space Y ( K ) . Thus, theelement problem for Maxwell’s equations is unisolvent.We prove that the global problem for Maxwell’s equations is unisolvent by showingthat ˆ E = ⃗ η = ˆ E in equation (16)and take the real part to get ∑ K ∈T h Im ( k ) (∣∣ ⃗ H ∣∣ K + ∣∣ ⃗ E ∣∣ K ) − Re ( τ )∣∣⃗ n × ( ˆ E − ⃗ E )∣∣ ∂K = . (19)This immediately shows that if condition (18b) holds, then the fields ⃗ E and ⃗ H are zero onΩ ⊂ R and the proof is finished. In the case of condition (18a), we have ⃗ n ×( ˆ E − ⃗ E ∣ ∂K ) = ⃗ K . Using equations (10), this yields [ ⃗∇×( ⃗∇× ⃗ E )]∣ K = k ⃗ E ∣ K , so Lemma 1 proves that the fields on element interiors are zero, which in turn impliesˆ E = ⃗ (cid:3) Note that even with Dirichlet boundary conditions and real k , the theorem asserts theexistence of a unique solution for the Helmholtz equation. However, the exact Helmholtzproblem (1) is well-known to be not uniquely solvable when k is set to one of an infinitesequence of real resonance values. The fact that the discrete system is uniquely solvableeven when the exact system is not, suggests the presence of artificial dissipation in HDGmethods. We will investigate this issue more thoroughly in the next section.However, we do not advocate relying on this discrete unisolvency near a resonancewhere the original boundary value problem is not uniquely solvable. The discrete ma-trix, although invertible, can be poorly conditioned near these resonances. Consider,for example, the Helmholtz equation on the unit square with Dirichlet boundary condi-tions. The first resonance occurs at k = π √
2. In Figure 2, we plot the condition number σ max / σ min of the condensed HDG matrix for a range of wavenumbers near the resonance k = π √
2, using a small fixed mesh of mesh size h = /
4, and a value of τ = p = p = TABILIZATION AND WAVENUMBER IN HDG 13 σ m a x / σ m i n τ = 1 τ = − i | (a) Degree p = σ m a x / σ m i n τ = 1 τ = − i | (b) Degree p = Figure 2.
Conditioning of the HDG matrix for the Helmholtz equationnear the first resonance k = π √ ≈ . τ = − ˆ ı that does not satisfy the conditions ofthe theorem produce much larger condition numbers, e.g., the condition numbers thatare orders of magnitude greater than 10 (off axis limits of Figure 2b) for k near theresonance were obtained for p = τ = − ˆ ı . To summarize the caveat, even thoughthe condition number is always bounded for values of τ that satisfy (18) , it may still bepractically infeasible to find the HDG solution near a resonance. Results of dispersion analysis for real wavenumbers
When the wavenumber k is complex , we have seen that it is important to choose thestabilization parameter τ such that (18b) holds. We have also seen that when k is real,the stability obtained by (18a) is so tenuous that it is of negligible practical value. Forreal wavenumbers, it is safer to rely on stability of the (un-discretized) boundary valueproblem, rather than the stability obtained by a choice of τ .The focus of this section is on real k and the Helmholtz equation (1). In this case,having already separated the issue of stability from the choice of τ , we are now freeto optimize the choice of τ for other goals. By means of a dispersion analysis, we nowproceed to show that some values of τ are better than others for minimizing discrepanciesin wavespeed. Since dispersion analyses are limited to the study of propagation ofplane waves (that solve the Helmholtz equation), we will not explicitly consider theMaxwell HDG system in this section. However, since we have written the Helmholtzand Maxwell system consistently with respect to the stabilization parameter (see thetransition from (3) to (9) via (7)), we anticipate our results for the 2D Helmholtz caseto be useful for the Maxwell case also. The dispersion relation in the one-dimensional case.
Consider the HDG method (2)in the lowest order ( p =
0) case in one dimension (1D) – after appropriately interpretingthe boundary terms in (2). We follow the techniques of [1] for performing a dispersion analysis. Using a basis on a segment of size h in this order u = , φ = , ˆ φ = , ˆ φ = , the HDG element matrix takes the form M = [ M M M M ] where M = [ ˆ ıkh − ˆ ıkh − τ ] M = [− + τ τ ] M = M t M = [− τ − τ ] . The Schur complement for the two endpoint basis functions { ˆ φ , ˆ φ } is then S = − ⎡⎢⎢⎢⎢⎢⎢⎢⎣ ıkh − τ ˆ ıkh + τ + τ − ıkh − τ ˆ ıkh + τ − ıkh − τ ˆ ıkh + τ ıkh − τ ˆ ıkh + τ + τ ⎤⎥⎥⎥⎥⎥⎥⎥⎦ . Applying this matrix on an infinite uniform grid (of elements of size h ), we obtain thestencil at an arbitrary point. If ˆ ψ j denotes the solution (trace) value at the j th point( j ∈ Z ), then the j th equation reads2 ˆ ψ j ( ıkh − τ ˆ ıkh + τ + τ ) + ( ˆ ψ j − + ˆ ψ j + ) (− ıkh − τ ˆ ıkh + τ ) = . In a dispersion analysis, we are interested in how this equation propagates plane waveson the infinite uniform grid. Hence, substituting ˆ ψ j = exp ( ˆ ık h jh ) , we get the followingdispersion relation for the unknown discrete wavenumber k h :cos ( k h h ) ( ıkh + τ ˆ ıkh + τ ) = ( ıkh − τ ˆ ıkh + τ + τ ) Simplifying, k h h = cos − ( − ( kh ) + ˆ ıkh ( τ + τ − ) ) . (20)This is the dispersion relation for the HDG method in the lowest order case in onedimension. Even when τ and k are real, the argument of the arccosine is not. HenceIm ( k h ) ≠ , (21)in general, indicating the presence of artificial dissipation in HDG methods . Note how-ever that if τ is purely imaginary and kh is sufficiently small, (20) implies that Im ( k h ) = kh (i.e., large number of elements per wavelength).As kh →
0, using the approximation cos − ( − x / ) ≈ x + x / + ⋯ valid for small x ,and simplifying (20), we obtain k h h − kh = − ( τ + ) ˆ ı τ ( kh ) + O (( kh ) ) . (22)Comparing this with the discrete dispersion relation of the standard finite elementmethod in one space dimension (see [1]), namely k h h − kh ≈ O (( kh ) ) , we find thatwavespeed discrepancies from the HDG method can be larger depending on the value of τ . In particular, we conclude that if we choose τ = ± ˆ ı , then the error k h h − kh in bothmethods are of the same order O (( kh ) ) . TABILIZATION AND WAVENUMBER IN HDG 15
Before concluding this discussion of the one-dimensional case, we note an alternateform of the dispersion relation suitable for comparison with later formulas. Using thehalf-angle formula, equation (20) can be rewritten as c = − ( ( kh ) ) ( τ ˆ ıkh ( τ + ) + τ ) , (23)where c = cos ( k h h / ) . Lowest order two-dimensional case.
In the 2D case, we use an infinite grid of squareelements of side length h . The HDG element matrix associated to the lowest order ( p = × C and C denotethe infinite set of stencil centers for the two types of stencils present in our case. Then,we get an infinite system of equations for the unknown solution (numerical trace) valuesˆ ψ , ⃗ p and ˆ ψ , ⃗ p at all ⃗ p ∈ C and ⃗ p ∈ C , respectively. We are interested in how thisinfinite system propagates plane wave solutions in every angle θ . Therefore, with theansatz ˆ ψ j, ⃗ p j = a j exp ( ˆ ı ⃗ κ h ⋅ ⃗ p j ) for constants a j ( j = ⃗ κ h = k h [ cos θ sin θ ] we proceed to find the relation between the discrete wavenumber k h and the exactwavenumber k .Substituting the ansatz into the infinite system of equations and simplifying, we obtaina 2 × F [ a a ] = F = [ khτ c c d ( τ + ˆ ıkh ) + khτ c d ( τ + ˆ ıkh ) + khτ c khτ c c ] and, for j = , c j = cos ( hk hj ) , d j = ı ( − c j ) − τ kh, k h = k h cos θ, k h = k h sin θ. (24)Hence the
2D dispersion relation relating k h to k in the HDG method isdet ( F ) = . (25)To formally compare this to the 1D dispersion relation, consider these two sufficientconditions for det ( F ) = ( kh ) τ c j + d j ( τ kh + ˆ ı ( k j h ) ) = , for j = , , (26) where k = k cos θ and k = k sin θ . (Indeed, multiplying (26) j by d j + ( j =
1) or d j − ( j =
2) and summing over j = ,
2, one obtains a multiple of det ( F ) .) The equationsin (26) can be simplified to c j = − ( k j h ) ı ( kh τ ( k j h ) + ( kh ) τ − ı kh τ ) , j = , , (27)which are relations that have a form similar to the 1D relation (23). Hence we useasymptotic expansions of arccosine for small kh , similar to the ones used in the 1D case,to obtain an expansion for k hj , for j = , k h = (( k h ) + ( k h ) ) / . (28)Simplifying the above-mentioned expansions for each term on the right hand side above,we obtain k h h − kh = ˆ ı ( cos ( θ ) + + τ ) τ ( kh ) + O (( kh ) ) (29)as kh →
0. Thus, we conclude that if we want dispersion errors to be O (( kh ) ) , then wemust choose τ = ± ı √ cos ( θ ) + , (30)a prescription that is not very useful in practice because it depends on the propagationangle θ . However, we can obtain a more practically useful condition by setting τ to bethe constant value that best approximates ± ˆ ı √ cos ( θ ) + ⩽ θ ⩽ π /
2, namely τ = ± ˆ ı √ . (31)These values of τ asymptotically minimize errors in discrete wavenumber over all anglesfor the lowest order 2D HDG method. Note that for any purely imaginary τ , (27) impliesthat k hj is real if kh is sufficiently small, soIm ( k h ) = , (32)thus eliminating artificial dissipation.We now report results of numerical computation of k h = k h ( θ ) by directly applyinga nonlinear solver to the 2D dispersion relation (25) (for a set of propagation angles θ ). The obtained values of the real part Re k h ( θ ) are plotted in Figure 3a, for a fewfixed values of τ . The discrepancy between the exact and discrete curves quantifies thedifference in the wave speeds for the computed and the exact wave. Next, analyzingthe computed k h ( θ ) for values of τ on a uniform grid in the complex plane, we foundthat the values of τ that minimize ∣ kh − k h ( θ ) h ∣ are purely imaginary. As shown inFigure 4, these τ -values approach the asymptotic values determined analytically in (30).A second validation of our analysis is performed by considering the maximum error overall θ for each value of τ and then determining the practically optimal value of τ . Theresults, given in Table 2, show that the optimal τ values do approach the analytically TABILIZATION AND WAVENUMBER IN HDG 17 π /4 π /20.80.850.90.9511.05 Angle θ W a v e s peed Exact wave speedDiscrete wave speed, τ =1Discrete wave speed, τ = − iDiscrete wave speed, τ =0.1+0.8i | (a) p = π /4 π /20.980.9911.01 Angle θ W a v e s peed Exact wave speedDiscrete wave speed, τ =1Discrete wave speed, τ = − iDiscrete wave speed, τ =0.1+0.9i | (b) p = Figure 3.
Numerical wave speed Re (⃗ k h ( θ )) as a function of θ for variouschoices of τ . Here, k = h = π / kh Optimal τ , Optimal τ ,Im ( τ ) > ( τ ) < π / . ı − . ıπ / . ı − . ıπ /
16 0 . ı − . ıπ /
32 0 . ı − . ıπ /
64 0 . ı − . ıπ /
128 0 . ı − . ıπ /
256 0 . ı − . ı Table 2.
Numerically found values of τ that minimize ∣ kh − k h ( θ ) h ∣ forall θ in the p = ± ˆ ı √ ≈ ± . ı . Further numerical results for the p = Higher order case.
To go beyond the p = k h = k h ( θ ) . The accompanyingdispersive, dissipative, and total errors are defined respectively by (cid:15) disp = max θ ∣ Re ( k h ( θ )) − k ∣ , (cid:15) dissip = max θ ∣ Im ( k h ( θ ))∣ ,(cid:15) total = max θ ∣ k h ( θ ) − k ∣ . (33)Again, we consider an infinite lattice of h × h square elements with the ansatz thatthe HDG degrees of freedom interpolate a plane wave traveling in the θ direction withwavenumber k h . The lowest order and next higher order HDG stencils are comparedin Figure 5. Note that the figure only shows the interactions of the degrees of freedom Angle θ ( I m τ ) (a) Im ( τ ) > kh= π /4kh= π /16kh= π /64(cos(4 θ )+3)/4 Angle θ ( I m τ ) (b) Im ( τ ) < Figure 4.
The values of τ that locally minimize ∣ kh − k h h ∣ are purelyimaginary. Here, ( Im ( τ )) is compared with asymptotic values (dashedlines). (a) (b) (c)(d) (e) (f) Figure 5.
Stencils corresponding to the shaded node types. ( A )–( B ):Two node types of the lowest order ( p =
0) method; ( C )–( F ): Four nodetypes of the first order ( p =
1) method.corresponding to the ˆ φ variable—the only degrees of freedom involved after eliminationof the ⃗ u and φ degrees of freedom via static condensation. The lowest order method hastwo node types (shown in Figures 5a–5b), while the first order method has four nodetypes (shown in Figures 5c–5f). For a method with S distinct node types, denote the TABILIZATION AND WAVENUMBER IN HDG 19 solution value at a node of the s th type, 1 ⩽ s ⩽ S , located at ⃗ lh ∈ R , by ψ s, ⃗ l . With ouransatz that these solution values interpolate a plane wave, we have ψ s, ⃗ l = a s e ˆ ı ⃗ k h ⋅⃗ lh , for some constants a s .Now, to develop notation to express each stencil’s equation, we fix a stencil within thelattice. Suppose that it corresponds to a node of the t th type, 1 ⩽ t ⩽ S , that is locatedat ⃗ h . For 1 ⩽ s ⩽ S , define J t,s = {⃗ l ∈ R ∶ a node of type s is located at (⃗ + ⃗ l ) h } and,for ⃗ l ∈ J t,s , denote the stencil coefficient of the node at location (⃗ + ⃗ l ) h by D t,s, ⃗ l . Thestencil coefficient is the linear combination of the condensed local matrix entries thatwould likewise appear in the global matrix of equation (17). Both it and the set J t,s aretranslation invariant, i.e., independent of ⃗ . Since plane waves are exact solutions to theHelmholtz equation with zero sources, the stencil’s equation is S ∑ s = ∑ l ∈ J t,s D t,s, ⃗ l ψ s, ⃗ +⃗ l = . Finally, we remove all dependence on ⃗ in this equation by dividing by e ˆ ı ⃗ k h ⋅⃗ h , so thereare S equations in total, with the t th equation given by S ∑ s = a s ∑ l ∈ J t,s D t,s, ⃗ l e ˆ ı ⃗ k h ⋅⃗ lh = . (34)Defining the S × S matrix F ( k h ) by [ F ( k h )] t,s = ∑ l ∈ J t,s D t,s, ⃗ l e ˆ ık h [ cos θ, sin θ ]⋅⃗ lh , we observe that non-trivial coefficients { a s } exist if and only if k h is such thatdet ( F ( k h )) = . (35)This is the equation that we solve to determine k h for a given θ for any order (cf. (25)).Results of the dispersion analysis are shown in Figures 3 and 6. These figures com-bine the results from previously discussed p = p = k = h = π /
4, i.e., 8 elements per wavelength. Figure 6shows the dispersive, dissipative, and total errors for various values of τ ∈ C . For boththe lowest order and first order cases, although the dispersive error is minimized at avalue of τ having nonzero real part, the total error is minimized at a purely imaginaryvalue of τ . This is attributed to the small dissipative errors for such τ . Specifically, thetotal error is minimized when τ = . ı in the p = τ found (both analytically and numerically) for p =
0. This value of τ reducesthe total wavenumber error by 90% in the p = τ = Re( τ ) I m ( τ ) − − − − (a) Dispersive error, p = Re( τ ) I m ( τ ) − − − − (b) Dispersive error, p = Re( τ ) I m ( τ ) − − − − (c) Dissipative error, p = Re( τ ) I m ( τ ) − − − − (d) Dissipative error, p = Re( τ ) I m ( τ ) − − − − (e) Total error, p = Re( τ ) I m ( τ ) − − − − (f) Total error, p = Figure 6.
Normalized dispersive error (cid:15) disp / (cid:15) , dissipative error (cid:15) dissip / (cid:15) , and total error (cid:15) total / (cid:15) for various τ ∈ C . Here, k = h = π /
4, and (cid:15) , (cid:15) and (cid:15) denote the errors when τ =
1, respec-tively.
TABILIZATION AND WAVENUMBER IN HDG 21
Comparison with dispersion relation for the Hybrid Raviart-Thomas method.
The HRT (Hybrid Raviart-Thomas) method is a classical mixed method [2, 3, 15] whichhas a similar stencil pattern, but uses different spaces. Namely, the HRT method forthe Helmholtz equation is defined by exactly the same equations as (2) but with thesechoices of spaces on square elements: V ( K ) = Q p + ,p ( K ) × Q p,p + ( K ) , W ( K ) = Q p ( K ) ,and M ( F ) = P p ( F ) . Here Q l,m ( K ) denotes the space of polynomials which are of degreeat most l in the first coordinate and of degree at most m in the second coordinate. Thegeneral method of dispersion analysis described in the previous subsection can be appliedfor the HRT method. We proceed to describe our new findings, which in the lowest ordercase includes an exact dispersion relation for the HRT method.In the p = k h for the HRT methodsatisfies the 2D dispersion relation ( c + c ) ( ( hk ) − ) + c c ( ( hk ) + ) + ( hk ) − = , (36)where c j , as defined in (24), depends on k hj , which in turn depends on k h . Similar to theHDG case, we now observe that the two equations ( ( hk j ) + ) c j + ( hk j ) − = , j = , , (37)are sufficient conditions for (36) to hold. Indeed, if l j is the left hand side above, then l ( c + ) + l ( c + ) equals the left hand side of (36). The equations of (37) canimmediately be solved: hk hj = − ( − ( hk j ) ( hk j ) + ) / Hence, using (28) and simplifying using the same type of asymptotic expansions as theones we previously used, we obtain k h h − kh = − cos ( θ ) + ( kh ) + O (( kh ) ) (38)as kh →
0. Comparing with (29), we find that in the lowest order case, the HRT methodhas an error in wavenumber that is asymptotically one order smaller than the HDGmethod for any propagation angle, irrespective of the value of τ .To conclude this discussion, we report the results from numerically solving the nonlin-ear solution (36) for k h ( θ ) for an equidistributed set of propagation angles θ . We havealso calculated the analogue of (36) for the p = h , these errors for boththe HDG and the HRT methods are graphed in Figure 7 for p = p =
1. We findthat the dispersive errors decrease at the same order for the HRT method and the HDGmethod with τ =
1. While (38) suggests that the dissipative errors for the HRT methodshould be of higher order, our numerical results found them to be zero (up to machineaccuracy). The dissipative errors also quickly fell to machine zero for the HDG methodwith the previously discussed “best” value of τ = ˆ ı √ /
2, as seen from Figure 7. −2 −1 −6 −4 −2 kh | k h − R e ( k h h ) | (a) Dispersive error, p = −2 −1 −15 −10 −5 kh | k h − R e ( k h h ) | (b) Dispersive error, p = −2 −1 −20 −10 kh | I m ( k h h ) | (c) Dissipative error, p = −2 −1 −20 −10 kh | I m ( k h h ) | (d) Dissipative error, p = −2 −1 −6 −4 −2 kh | k h − k h h | (e) Total error, p = −2 −1 −15 −10 −5 kh | k h − k h h | (f) Total error, p = Order (kh) referenceOrder (kh) referenceOrder (kh) referenceOrder (kh) referenceHDG with τ = 1HDG with τ = 0.866 iHRT Figure 7.
Convergence rates as kh → TABILIZATION AND WAVENUMBER IN HDG 23
Conclusions
These are the findings in this paper:(1) There are values of stabilization parameters τ that will cause the HDG methodto fail in time-harmonic electromagnetic and acoustic simulations using complexwavenumbers. (See equation (5) et seq.)(2) If the wavenumber k is complex, then choosing τ so that Re ( τ ) Im ( k ) ⩽ k is real, then even when the exact wave problem is not well-posed (such as at a resonance), the HDG method remains uniquely solvable whenRe ( τ ) ≠
0. However, in such cases, we found the discrete stability to be tenuous.(See Figure 2 and accompanying discussion.)(4) For real wavenumbers k , we found that the HDG method introduces smallamounts of artificial dissipation (see equation (21)) in general. However, when τ is purely imaginary and kh is sufficiently small, artificial dissipation is eliminated(see equation (32)). In 1D, the optimal values of τ that asymptotically minimizethe total error in the wavenumber (that quantifies dissipative and dispersive er-rors together) are τ = ± ˆ ı (see equation (22)).(5) In 2D, for real wavenumbers k , the best values of τ are dependent on the propa-gation angle. Overall, values of τ that asymptotically minimize the error in thediscrete wavenumber (considering all angles) is τ = ± ˆ ı √ / τ = ˆ ı √ /
2, dissipative errorsdominate when τ = References [1]
M. Ainsworth , Discrete dispersion relation for hp -version finite element approximation at highwave number , SIAM J. Numer. Anal., 42 (2004), pp. 553–575 (electronic).[2] D. N. Arnold and F. Brezzi , Mixed and nonconforming finite element methods: implementa-tion, postprocessing and error estimates , RAIRO Mod´el. Math. Anal. Num´er., 19 (1985), pp. 7–32.[3]
B. Cockburn and J. Gopalakrishnan , A characterization of hybridized mixed methods for theDirichlet problem , SIAM J. Numer. Anal., 42 (2004), pp. 283–301.[4]
B. Cockburn, J. Gopalakrishnan, and R. Lazarov , Unified hybridization of discontinuousGalerkin, mixed, and continuous Galerkin methods for second order elliptic problems , SIAM Journalon Numerical Analysis, 47 (2009), pp. 1319–1365.[5]
B. Cockburn, J. Gopalakrishnan, and F.-J. Sayas , A projection-based error analysis ofHDG methods , Math. Comp., 79 (2010), pp. 1351–1367.[6]
J. Cui and W. Zhang , An analysis of HDG methods for the Helmholtz equation , IMA J. Numer.Anal., 34 (2014), pp. 279–295.[7]
A. Deraemaeker, I. Babuˇska, and P. Bouillard , Dispersion and pollution of the FEMsolution for the Helmholtz equation in one, two and three dimensions , International Journal forNumerical Methods in Engineering, 46 (1999), pp. 471–499. [8]
G. Giorgiani, S. Fern´andez-M´endez, and A. Huerta , Hybridizable discontinuous Galerkinp-adaptivity for wave propagation problems , International Journal for Numerical Methods in Fluids,72 (2013), pp. 1244–1262.[9]
J. Gopalakrishnan, I. Muga, and N. Olivares , Dispersive and dissipative errors in the DPGmethod with scaled norms for the Helmholtz equation , SIAM J. Sci. Comput., 36 (2014), pp. A20–A39.[10]
R. Griesmaier and P. Monk , Error analysis for a hybridizable discontinuous Galerkin methodfor the Helmholtz equation , J. Sci. Comput., 49 (2011), pp. 291–310.[11]
A. Huerta, X. Roca, A. Aleksandar, and J. Peraire , Are high-order and hybridizablediscontinuous Galerkin methods competitive? , in Oberwolfach Reports, vol. 9 of Abstracts fromthe workshop held February 12–18, 2012, organized by Olivier Allix, Carsten Carstensen, J¨orgSchr¨oder and Peter Wriggers, Oberwolfach, Blackforest, Germany, 2012, pp. 485–487.[12]
L. Li, S. Lanteri, and R. Perrussel , Numerical investigation of a high order hybridizablediscontinuous Galerkin method for 2d time-harmonic Maxwell’s equations , COMPEL, 32 (2013),pp. 1112–1138.[13] ,
A hybridizable discontinuous Galerkin method combined to a Schwarz algorithm for thesolution of 3d time-harmonic Maxwell’s equation , J. Comput. Phys., 256 (2014), pp. 563–581.[14]
N. Nguyen, J. Peraire, and B. Cockburn , Hybridizable discontinuous Galerkin methods forthe time-harmonic Maxwell’s equations , Journal of Computational Physics, 230 (2011), pp. 7151–7175.[15]
P.-A. Raviart and J. M. Thomas , A mixed finite element method for 2nd order elliptic problems ,in Mathematical aspects of finite element methods (Proc. Conf., Consiglio Naz. delle Ricerche(C.N.R.), Rome, 1975), Springer, Berlin, 1977, pp. 292–315. Lecture Notes in Math., Vol. 606.
Portland State University, PO Box 751, Portland, OR 97207-0751, USA
E-mail address : [email protected] INRIA Sophia Antipolis M´editerran´ee, 2004 Route des Lucioles, BP 93, 06902 SophiaAntipolis Cedex, France
E-mail address : [email protected] Portland State University, PO Box 751, Portland, OR 97207-0751, USA
E-mail address : [email protected] LAPLACE (LAboratoire PLAsma et Conversion dEnergie), Universit´e de Toulouse,CNRS/INPT/UPS, Toulouse,France
E-mail address ::