A comparison of high-order polynomial and wave-based methods for Helmholtz problems

TL;DR: The high-order polynomial method (p-FEM with Lobatto polynomials) and the wave-based discontinuous Galerkin method are compared for two-dimensional Helmholtz problems, indicating that the differences in performance, accuracy and conditioning are more nuanced than generally assumed.

About: This article is published in Journal of Computational Physics.The article was published on 2016-09-15 and is currently open access. It has received 44 citations till now. The article focuses on the topics: Polynomial & Helmholtz equation.

Summary

1. Introduction

• Numerical dispersion is the main source of inaccuracy in the conventional finite element method when solving Helmholtz problems at high frequencies.
• It is generally assumed that a plane-wave basis naturally provides better accuracy compared to polynomials since the latter bears no relation with the governing equations.
• A number of wave-based methods rely on discontinuous formulations where the solution is approximated with plane waves in each element and the continuity between elements is imposed weakly.
• Numerical results in [19, 26] indicate that it is more efficient than the UWVF and least-square methods.

2. Presentation of the methods

• On the remainder of the boundary EQUATION 2.1.
• This section describes the formulation of the high-order, polynomial finite-element method (p-FEM) for solving the Helmholtz equation.

2.1.2. High-order shape functions

• The approximate solution u h and the associated test function v h are constructed using the classical H 1 -conforming hierarchical polynomial shape functions.
• This reference element is equipped with the following affine coordinates: EQUATION ).
• The construction of the higher order shape functions rely on the Lobatto polynomials EQUATION with L q the Legendre polynomial of order q ≥.
• These functions verify the following orthogonality property: 1 EQUATION which ensures an optimal conditioning of the stiffness matrix for Helmholtz problems [38] .
• For the evaluation of the integrals in the p-FEM formulation (4), Gauss-Legendre quadratures are used [6] .

2.1.3. Static condensation

• The bubble functions defined in (11) vanish on the element boundaries and are therefore not coupled to the neighbouring elements.
• These internal degrees of freedom can be eliminated at the element level prior to assembly and do not appear in the global system matrix.
• This static condensation procedure does not affect the final solution and is a very effective technique to reduce the size and improve the conditioning of the global system matrices [6, 7] .
• The benefits of condensation tends to become more significant as the polynomial order p increases because the number of bubble functions scales like p 2 whereas the number of edge functions scales like p.
• Note however that linear and quadratic triangular elements do not involve bubble functions and cannot benefit from static condensation.

2.1.4. Error estimates

• The first term in this expression represents the interpolation error and has a dependency of order p with kh.
• This error can be controlled by using a sufficient number of elements per wavelength, given by 2π/(kh).
• The second term is associated with the dispersion error and the pollution effect.
• From this error estimate, the advantages of resorting to higher order shape functions appears clearly for smooth solutions.
• Exponential convergence rates are obtained with p-refinement, while convergence with h-refinement is only algebraic.

2.2. Wave-based DGM

• This section describes the formulation of the wave-based DGM [18] for solving the Helmholtz equation.
• Only the key aspects of the method that are required for the discussion of the results are described here.
• The reader is referred to [19] for a more detailed presentation.

2.2.2. Variational Formulation

• After integration by parts in each element the variational formulation for equation (13) reads: EQUATION where T denotes the conjugate transpose.
• In the wave-based DGM, the solution in each element is assumed to satisfy the homogeneous equation: EQUATION.

2.2.3. Numerical flux

• The different categories of numerical fluxes available will not be reviewed here and the reader is referred to [42] .
• In the current formulation, the upwind flux-vector splitting method is used.
• The flux on the edge between two elements is expressed in terms of characteristic waves in the direction normal to the edge [43] .
• The expressions for Λ and W are available in [19] and are not repeated here.

2.2.4. Boundary condition

• In the present paper, ghost cells and Dirichlet boundary conditions are used.
• This technique follows naturally from the numerical flux presented above and the idea is to treat the outer boundary of the computational domain as an internal edge between two elements.
• It should be noted that this approach is equivalent to the Robin boundary condition (2) used for the p-FEM.
• For the Dirichlet boundary condition (3), the method outlined in [44] is used to modify the variational formulation accordingly.
• Imposing the pressure on the boundary amounts to specifying a relation between the outgoing and incoming characteristic waves.

2.2.5. Interpolation basis

• With Trefftz methods the solution is approximated using local, canonical solutions of the governing equations.
• For simplicity, in this study the same number N w of plane waves is used in every element and the wave directions are evenly distributed on the interval [0, 2π[.
• While not shown here, results for even numbers of plane waves have also been calculated and the performance and conditioning of the wave-based DGM was indeed found to be slightly below those obtained with odd values of N w .
• By construction the approximate solution will directly include these properties.

2.2.6. Error estimates

• Theoretical error estimates for the h-version of the UWVF are given in [51] , for the Helmholtz equation.
• However numerical tests carried out for different frequencies and number of plane waves indicate that this expression under estimates the actual convergence rates.
• Similar error estimates are also derived for h-refinement by Cessenat et al. in [14] , although for a slightly different error norm.
• In a series of papers, Hiptmair et al. [20, 22, 23] discussed the behavior of the p and the hp-version of the wavebased DGM.
• The relative dispersion error is shown to behave asymptotically like (kh) (N w −1) for large values of kh, hence exhibiting an exponential decrease as the number of plane waves is increased.

3. Description of the test cases

• The test cases used to assess the p-FEM and the wave-based DGM are introduced.
• The first test case is a simple plane wave propagating on a square domain (see Fig. 1(a) ).
• All the benchmark problems described above involve smooth solutions, but the performance of high-order schemes is known to deteriorate significantly when dealing with non-smooth solutions, including the p-FEM [53] and the wavebased DGM [36] .
• Therefore, the fourth test case involves a corner solution.
• For all the other boundaries the Robin condition ( 2) is used for p-FEM and the ghost cells (23) are used for the wave-based DGM.

4. Measures of accuracy and costs

• For all the test cases described above a number of parameters can be varied: the angular frequency ω, the element size h and the approximation order (that is the polynomial order p for the p-FEM and the number of plane waves N w for the wave-based DGM).
• Unless explicitly stated the mesh resolution will be measured using D λ instead of the element size h.
• Furthermore, the matrices arising from a high-order discretisation typically exhibit higher condition numbers than that of the linear FEM, which renders the problem even more challenging (high-order bases designed to address this issue have been proposed [55, 56] ).
• Its effects on computational cost and conditioning will also be addressed.

5. Comparison of the interpolation properties

• Before trying to solve the test cases previously defined using the two methods, a first analysis is performed to assess the interpolation properties of their respective approximation bases.
• Independently of the numerical formulations, this section discusses the ability of the continuous polynomial basis on the one hand, and the discontinuous plane-wave basis on the other hand, to approximate the solutions of the four different test cases.
• These two integrals are also evaluated using high-order Gauss-Legendre quadrature rules, with orders 2p + 10 and N w +10 for the p-FEM and the wave-based DGM methods, respectively.
• The relative error on the optimal interpolation is denoted by EQUATION.
• It is worth noting that the results for the wave-based DGM presented in this section also apply to a number of other wave-based methods: UWVF, least-squares method and wave-based DGM with Lagrange multipliers, since all of these methods rely on the same plane-wave approximation basis defined in Section 2.2.5.

5.1. Anisotropy

• The anisotropy of the two different approximating bases is best shown using the test case of the single plane wave by varying the incident wave direction θ 0 from 0 to 2π.
• The behaviour of the two bases with respect to the propagation angle θ 0 is very different.
• In contrast, the plane-wave basis is very sensitive to the orientation of the incident plane wave.
• When it is aligned with one of the waves in the basis, the error drops to zero since the exact solution belongs to the interpolation basis [18] .
• This strong anisotropy of the wave-based DGM also justifies the introduction of the spinning wave test case ) where all wave directions are equally represented.

5.2. Convergence properties

• The interpolation properties of the two approximation bases are now examined with respect to order and mesh resolution.
• For the single plane wave and the propagating spinning wave test cases, the expected behaviour is observed, with an algebraic convergence with element size h and an exponential convergence with p or N w .
• This is due to the poor conditioning of the system solved to obtain the best interpolation.
• For the wave-based DGM the numerical error is found to be concentrated close to the cylinder which is consistent with the evanescent part of the solution being less accurately resolved by the plane-wave basis, as shown in figure 4(c) .
• This approach obviously requires a detailed knowledge of the solution behaviour before the numerical method is even formulated.

6. Performance of the numerical methods

• The efficiency of the two numerical methods is compared.
• The convergence with respect to p-and h-refinements is firstly considered, followed by the study of the conditioning properties of the underlying system of equations.
• Finally, the computational cost at fixed accuracy is assessed.

6.1. Convergence

• For the plane wave and the spinning wave test cases, these convergence results follow, at least qualitatively, the results in Figure 4 for the best interpolation error.
• High-order methods allow to mitigate the pollution effect and this is indeed observed in Figure 6 .
• It can be noted that, overall, p-FEM exhibits slightly better optimality coefficients.
• This indicates that not only is the corner solution an issue in terms of interpolation but the expected convergence properties of high-order methods are also lost.
• These results illustrate that h-refinement has a significantly different impact on the behaviour of the two methods.

6.2. Conditioning

• In Figure 9 the condition numbers are plotted against D λ for the propagating spinning wave test case.
• The impact of static condensation on the conditioning is very significant.
• Without static condensation, the conditioning of the p-FEM model deteriorates very rapidly when the polynomial order is increased.
• In particular, for a large number of plane waves, the condition number grows very steeply with the mesh resolution.
• Hence, for reasonable levels of accuracy, say for instance 1%, the conditioning of the wave-based method is not problematic, and is in fact equivalent to the p-FEM with static condensation.

6.3. Comparison of computational costs

• Given the large number of parameters involved, it is convenient to compare the two methods for a fixed accuracy and frequency.
• One can assess the requirements needed for the two methods to solve for a particular problem with a given error threshold.
• In a second stage, the authors will discuss to which extent the conclusions drawn are applicable to other frequencies and other levels of accuracy.
• For a given test case and for every order p or number of plane waves N w , the mesh size h is varied until the numerical error reaches the target precision, then the cost and conditioning of this model are recorded.
• The test case of the spinning wave is used here but the conclusions are equally applicable to the plane-wave test case.

6.3.1. Fixed accuracy and frequency

• Figure 10 shows the link between the mesh size and the factorisation time.
• Also shown in Figure 10 is the reduction in factorisation time obtained from static condensation (the mesh resolution remains the same with and without static condensation).
• For the propagating spinning wave, see Figure 10 (a), increasing N w for the wave-based method leads to a consistent reduction in the memory required to solve the problem but the conditioning does not vary significantly.
• For p-FEM with static condensation, increasing the polynomial order also reduces the factorisation memory and induces a moderate increase in conditioning, although the condition number remains comparable to that of the wave-based method.
• This difference is further accentuated when considering the evanescent spinning wave.

6.3.2. Influence of the target accuracy

• The test case of the propagating spinning wave is used, but the results for the plane wave and the evanescent spinning wave test cases are similar.
• For the factorisation memory, the two methods follow the same trends: the memory increases when the target error level is reduced and when D λ increases.
• The fact that the wave-based method requires slightly more memory than p-FEM is observed at all error levels, confirming that this observation is not specific to the 1% target accuracy used above.
• For the wave-based method, the condition number increases very rapidly to reach levels that are problematic.
• Again, this is consistent with the results in Figure 9 which shows that it is only for very accurate solutions that the conditioning of the wave-based method becomes a real issue.

6.3.3. Influence of the frequency

• Finally, it is important to assess whether the trends described above are also observed for other frequencies.
• To this end, Figure 14 shows the factorisation memory and the conditioning for Helmholtz numbers ranging from kL = 25 to 200 for the propagative spinning wave test case (the azimuthal order m is adjusted at each frequency to ensure the wave is propagating).
• As the polynomial order or the number of plane waves is increased, the two methods are able to achieve 1% error with less memory.
• For every frequency, the wave-based method requires a larger amount of memory than p-FEM, although the difference is relatively small for the lowest frequency (kL = 25).
• The evolution of the condition number with respect to the frequency is more complex and does not follow a simple trend.

7. Practical considerations

• So far, the discussion has focused on the intrinsic computational performance of the numerical schemes.
• Both methods lose p-convergence when dealing with singular solutions.
• The presence of source terms in the domain also requires special techniques to address the fact that the plane waves are solutions of the homogeneous equations.
• While this is acceptable for applications to academic test cases, it is difficult to envisage for real-world applications where robust, fully automated solvers are required.
• By contrast, with p-FEM the use of hierarchic shape functions allows to store and reuse the element matrices for other frequencies [10] .

8. Conclusions

• A detailed comparison of a high-order polynomial method (p-FEM) and a wave-based method (DGM) has been presented in terms of interpolation properties, performance and conditioning.
• The main objective was to put in perspective the relative merits of these two categories of discretisations.
• The general argument for the use of wave-based methods is that plane waves are canonical solutions of the underlying equations and are therefore better suited than polynomials to construct an approximation basis.
• This observation was found to hold for a range of frequencies and error levels.
• It should be noted that the use of static condensation is crucial to obtain good performance and conditioning from p-FEM.

Figures

Figure 5: Relative L2 error (%) against Dλ at kL = 50; red dashed lines: wave-based DGM; black solid lines: p-FEM with condensation. For each curve, Dλ is increased by refining the mesh while p or Nw is fixed.

Figure 12: Number of non-zero entries against the number of degree of freedom for a numerical error of 1% at kL = 50. Red dashed lines: wave-based DGM; Blue dotted lines: p-FEM without condensation; Black solid lines: p-FEM with condensation.

Figure 13: Propagating spinning wave test case (kL = 50, m = 10); factorisation memory (left) and condition number (right) against Dλ to achieve different levels of accuracy; red dashed lines: wave-based DGM, solid black lines: p-FEM with static condensation.

Figure 14: Propagating spinning wave problem. Factorisation memory (left) and condition number (right) against Dλ to achieve 1% of accuracy at different frequencies; red dashed lines: wave-based DGM, solid black lines: p-FEM with static condensation.

Figure 9: Propagating spinning wave case with kL = 50 and m = 10. Left: Condition number with respect to Dλ (For each curve, Dλ is increased by refining the mesh while p or Nw is fixed). Right: relative L2 error with respect to the condition number. Black solid lines: p-FEM with condensation, blue dotted lines: p-FEM without condensation, red dashed lines: wave-based DGM.

Figure 4: Relative best L2 interpolation error (%) against Dλ at kL = 50; red dashed lines: wave-based method; black solid lines: p-FEM (with static condensation). For each curve, Dλ is increased by refining the mesh while p or Nw is fixed.

Figure 10: Factorisation time (s) against the element size required to obtain a numerical error of 1% at kL = 50. Red dashed lines: wave-based DGM; Blue dotted lines: p-FEM without condensation; Black solid lines: p-FEM with condensation.

Figure 8: Relative L2 error (left) and condition number (right) against Cr at kL = 50 with hmax/b = 0.17 for the non-smooth problem with locally refined meshes; red dashed lines: wave-based DGM; black solid lines: p-FEM with condensation.

Figure 7: Examples of unstructured non-uniform meshes used for the non-smooth test case with hmax/b = 0.17.

Figure 3: Relative L2 error associated with the best interpolation plotted against the direction of the incident plane wave with h/a = 0.1, kL = 20. Black solid lines: p-FEM, red dashed lines: wave-based DGM. Numbers indicates the polynomial order or the number of plane waves.

Figure 1: Examples of analytical solutions (real part) of the four benchmark problems (kL = 50).

Figure 2: Examples of unstructured meshes used for the test cases.

Figure 11: Factorisation memory against condition number for a numerical error of 1% at kL = 50. Red dashed lines: wave-based DGM; Blue dotted lines: p-FEM without condensation; Black solid lines: p-FEM with condensation.

Figure 6: Optimality coefficient against Dλ for the propagating spinning wave with kL = 50; red lines: wave-based DGM; black lines: p-FEM with condensation. For each curve, Dλ is increased by refining the mesh while p or Nw is fixed.

Frequently asked questions

Q1. What are the contributions in "A comparison of high-order polynomial and wave-based methods for helmholtz problems" ?

The application of computational modelling to wave propagation problems is hindered by the dispersion error introduced by the discretisation. In this paper a high-order polynomial method ( p-FEM with Lobatto polynomials ) and the wave-based discontinuous Galerkin method are compared for two-dimensional Helmholtz problems. 

Q2. What is the main source of inaccuracy in the conventional finite element method when solving?

Numerical dispersion is the main source of inaccuracy in the conventional finite element method when solving Helmholtz problems at high frequencies. 

Q3. What is the general argument for the use of wave-based methods?

The general argument for the use of wave-based methods is that plane waves are canonical solutions of the underlying equations and are therefore better suited than polynomials to construct an approximation basis. 

Q4. What is the way to improve the FEM and wave-based DGM solutions?

The use of local h-refinement is a valid option to improve both the FEM and the wave-based DGM solutions and will be used in the next section. 

Q5. How many tests have been performed without significant change in the results?

2. A large number of tests have also been performed using quadrature orders up to 2p + 10, without inducing any significant change in the results, thereby indicating that the resulting integration error can be considered negligible. 

Q6. What is the central argument for using a basis of plane waves?

The central argument for using a basis of plane waves is that it follows directly from the governing equations and thus it allows embedding key features of the underlying physics into the numerical method. 

Q7. What is the first step in the study of the underlying system of equations?

The convergence with respect to p- and h-refinements is firstly considered, followed by the study of the conditioning properties of the underlying system of equations. 

Q8. What is the continuity of the normal flux across the edge between the elements e and e?

The continuity of the normal flux across the edge between the elements e and e′ is directly enforced by writing:Feue = −Fe′ue′ = fe,e′ (ue,ue′ ) , (17)where fe,e′ is the so-called numerical flux which is discussed below. 

Q9. What is the generalized discretization error for the L2-norm?

It is shown that for smooth solutions, the global discretization error measured in the L2-norm behaves asymptotically like h(Nw−1)/2−1. 

Q10. How many different benchmark problems are used to assess the performance of high-order methods?

The performance of high-order methods are known to be problem dependent [12], therefore a total of four different benchmark problems are used to assess the ability of these methods to tackle a large variety of problems. 

Q11. What is the numerical method for constant coefficient matrices?

The numerical method is described here for constant coefficient matrices but it could be extended to the case where the matrices are functions of x and y [18] by approximating them as piecewise constant functions. 

Q12. What is the approximate solution uh and the associated test function vh?

The approximate solution uh and the associated test function vh are constructed using the classical H1-conforming hierarchical polynomial shape functions.