Adaptive PU-BEM for Helmholtz problems

1. Adaptive PU-BEM for Helmholtz problems Jon Trevelyan & Graham Coates School of Engineering, Durham University

2. The motivation Short wave scattering Limitation of polynomial basis elements (FEM/BEM) – 10 DOF per wavelength Partition of Unity Method (PUM) is just one approach that aims to overcome the limitation. The PUM is a general class of methods involving enrichment of the approximation space with sets of analytical functions known to populate the solution space. Mafelap Brunel University, June 2009

3. Partition of Unity Method Melenk & Babuška, Comp. Meth. Appl. Mech. Engg., 1996 Application to waves – PU-FEM Laghrouche & Bettess, J. Comput. Acoust., 2000. Ortiz & Sanchez, Int. J. Num. Meth. Eng., 2000. Farhat et al., various, discontinuous enrichment method Gamallo & Astley, Int. J. Num. Meth. Eng., 2006. Strouboulis et al, Comp Meth Appl Mech Eng., 2006. Ultraweak variational formulation Cessenat & Despres, SIAM J. Numer. Anal., 1998. Mafelap Brunel University, June 2009

4. Partition of Unity Method in integral equations Enrichment using plane wave(s) in integral equations de la Bourdonnaye, CRAS, 1994. “microlocal discretisation” Abboud et al., SIAM Wave propagation conf., 1995. Perrey-Debain et al., J. Sound & Vib., 2003. Bruno et al., Phil Trans Royal Society A, 2003. Langdon & Chandler-Wilde, SIAM J. Numer. Anal., 2006. etc. Mafelap Brunel University, June 2009

5. Incident plane wave The scattering problem In the frequency domain we solve the Helmholtz equation is the potential we seek and the wave number, k, is given by Mafelap Brunel University, June 2009

6. The scattering problem In the frequency domain we solve the Helmholtz equation is the potential we seek and the wave number is given by For simplicity consider Neumann problem: Mafelap Brunel University, June 2009

7. Collocation boundary integral equation Boundary integral equation collocating at x0 on smooth boundary  where the double layer potential operator is and the Green’s function (for 2D problems) is Mafelap Brunel University, June 2009

8. PU-BEM enriched basis The PUM multiple plane wave expansion for potential on element Polynomial shape functions Plane waves Amplitudes Reformulates the problem so that the unknowns become the amplitudes. Mafelap Brunel University, June 2009

9. PU-BEM Substituting the basis into the boundary integral equation and collocating at a sufficient number of points x0 gives a system of equations that may be solved for the amplitudes. To overcome the BIE non-uniqueness problem we use the CHIEF method of Schenck. Mafelap Brunel University, June 2009

10. Collocation method Each element will now have many degrees of freedom. degrees of freedom > number of nodes So we generate an auxiliary set of equations by collocating at a set of points uniformly distributed over the element. Mafelap Brunel University, June 2009

11. PU-BEM efficiency Define measure of efficiency  = DOF per wavelength Conventional polynomial basis requires   10 PU-BEM plane wave basis requires   2.5 Computational burden shifts from solver to integration and assembly. Important to optimise plane wave basis to minimise number of evaluations of highly oscillatory integrals. Mafelap Brunel University, June 2009

12. The number of plane waves at each node is increased as suggested by an error indicator Adaptive scheme We develop an adaptive scheme, of p-adaptive character, by allowing the number of plane waves at each node to vary. With the PUM basis our BIE becomes: Mafelap Brunel University, June 2009

13. Adaptive scheme BIE Define a residual error indicator Mafelap Brunel University, June 2009

14. Clearly where x1 coincides with one of the original collocation points, R will be very small: 0.015 R 0.01 0.005 0 x -1 -0.5 0 0.5 1 Variation in R over an element containing 13 collocation points Adaptive scheme Implementation of residual error indicator Evaluate at any point Mafelap Brunel University, June 2009

15. 0.015 R 0.01 0.005 0 x -1 -0.5 0 0.5 1 Adaptive scheme R has useful local and global properties. Global error indicator gives us a stopping criterion Define a 1-norm in terms of the value of R over a scatterer of perimeter P : This is expensive to compute, so we approximate using Mafelap Brunel University, June 2009

16. Adaptive scheme R has useful local and global properties. Global error indicator gives us a stopping criterion Numerical tests suggest a stopping criterion of gives a converged solution of engineering accuracy (1% relative error norm in L2) Mafelap Brunel University, June 2009

17. Adaptive scheme Local properties of R are used to decide where to add waves. To verify this, plot norms of both R and L2 error, evaluated over each element: 1.2 Norm of R and L2 error 1 normalised Rnorm 0.8 normalised L2 0.6 0.4 0.2 0 0 4 8 12 16 20 24 Element Mafelap Brunel University, June 2009

18. Adaptive scheme We evaluate R only at points close to end nodes 0.015 R 0.01 0.005 0 x -1 -0.5 0 0.5 1 Numerical tests suggest addition of a plane wave at nodes adjacent to points where R > 0.015 Mafelap Brunel University, June 2009

19. Plane wave addition At each adaptive step, we can introduce an extra plane wave at selected nodes. Mafelap Brunel University, June 2009

20. uniform uniform + 1 uniform + 2 Plane wave addition Check accuracy of PU-BEM for non-uniformly distributed wave directions. 12 elements Mafelap Brunel University, June 2009

21. Plane wave addition Check accuracy of PU-BEM for non-uniformly distributed wave directions. 0 uniform -1 uniform + 1 uniform + 2 -2 -3 -4 -5 -6 0 0.5 1 1.5 2 Mafelap Brunel University, June 2009

22. Adaptive scheme The extra DOF require us to add extra rows and columns A x = b New terms in blue areas to be found from boundary integrals Big savings in numerically intensive oscillatory integrals Mafelap Brunel University, June 2009

23. Illustration Scattering of a plane wave by sound-hard cylinder Initial model: 24 boundary elements 48 nodes 6 waves/node 288 DOF  = 2.29 a= 10 + l= 0.5 ka= 125 Mafelap Brunel University, June 2009

24. 0.07 R 0.06 0.05 0.04 0.03 0.02 0.01 0  0 2 4 6 Illustration  first run : 288 DOF Error indicator suggests a wave to be added at 21 nodes Mafelap Brunel University, June 2009

25. 0.07 R 0.06 0.05 0.04 0.03 0.02 0.01 0  0 2 4 6 Illustration  R 0.06 0.05 0.04 0.03 0.02 0.01 0  0 2 4 6 second run : 309 DOF first run : 288 DOF Error indicator suggests a wave to be added at 21 nodes Mafelap Brunel University, June 2009

26.  Illustration Converged solution Numerical Analytical  Mafelap Brunel University, June 2009

27. Second run Matrix 357  309  110313 matrix terms re-use 96768 13545 new terms Illustration Initial run Matrix 336  288  96768 matrix terms The other alternative without adaptivity… … run new analysis with 7 waves/node requiring 129024 terms Mafelap Brunel University, June 2009

28. Conclusions Adaptive scheme for PU-BEM in wave scattering Residual error indicator R Add waves near nodes having high R Norm of R gives stopping criterion Illustration shows large savings in required number of oscillatory integrals Mafelap Brunel University, June 2009