Modulated plane wave methods for Helmholtz problems in heterogeneous media Timo Betcke t.betcke@ucl.ac.uk Department of Mathematics University College London November 21, 2011 Joint work with Alex Barnett (Dartmouth), Joel Phillips (UCL) Supported by EPSRC Grant EP/H004009/1
Helmholtz Problems • Trefftz methods for homogeneous media � ω � 2 • Polynomially modulated plane ∆ u + u = 0 c ( x ) waves in inhomogeneous media ω k = c • Combination with high-frequency methods
From low to high k • Low k : Standard boundary and finite elements. Almost anything works. • k → ∞ : Raytracing, geometric theory of diffraction, fast eikonal solvers [Engquist/Runborg ’03]. • Midfrequency regime: • Speeding up the linear algebra: FMM [Cheng et. al. ’06, Engquist/Ying ’07], Novel preconditioners Engquist/Ying ’11], • Nonstandard basis functions: PUFEM [Babuska/Melenk ’97], UWVF [Cessenat/Despres ’98, Huttunen/Monk/Kaipio ’02], MPS/MFS [Barnett/B. ’10], Plane Wave DG Methods [Gittelson et. al ’09]
The Trefftz framework Elliptic PDE: Lu = 0 ( L := − ∆ − k 2 ). � Ω := K i : Decomposition of domain Ω into elements Ω i . i n α ( i ) � Local approximation spaces V i := { j φ j : L φ j = 0 in K i } . j = 1 • Local vs. global methods • Choices of basis functions • Assembly of discrete problem
Trefftz basis functions Plane Waves φ j = e ikd j · x , � d j � = 1 � cos θ j � , θ j = 2 π ( j − 1 ) d j = sin θ j N Fourier-Bessel functions φ j ( r , θ ) = J j ( kr ) e ij θ j = − n , . . . , n
Trefftz basis functions . . . Fundamental Solutions φ j ( x ) = H ( 1 ) 0 ( k | x − y j | ) y j : Source points • Plane waves can be chosen according to raytracing directions • Bessel functions have similar approximation properties as polynomials • Fundamental solutions good for approximating far-field and circular wavefronts • Differences in numerical stability between these bases [B. ’08, Barnett/B. ’08]
Approximation properties � 1 Vekua operator (2d version): V j [ φ ]( x ) = φ ( x ) + M j ( x , t ) φ ( tx ) dt 0 √ k | x | M 1 ( x , t ) = − √ 1 − t J 1 ( k | x | 1 − t ) 2 k | x | � M 2 ( x , t ) = − i √ √ J 1 ( ik | x | t ( 1 − t )) 2 t 1 − t Properties: 1 V 1 [ V 2 [ φ ]] = V 2 [ V 1 [ φ ]] = φ . 2 If φ is harmonic, then ∆ V 1 [ φ ] + k 2 V 1 [ φ ] = 0 . 3 If ∆ u + k 2 u = 0, then ∆ V 2 [ u ] = 0 [Eisenstat ’74, Still ’89, Melenk ’99, B. ’07, Moiola/Hiptmair/Perugia ’11]
Approximation properties . . . Generalized harmonic polynomials n � a j r | j | e ij φ P ( r , θ ) = j = − n n � | l | � 2 � e ij θ J | j | ( kr ) V 1 [ P ]( r , θ ) = a j | j | ! k j = − n n n � α j e id j · x � ≤ � V 1 �� φ − P � + � V 1 [ P ] − � α j e id j · x � � u − j = 1 j = 1 • Approximation error between plane waves and Fourier-Bessel fct • Polynomial approximation of harmonic functions Exponential estimates via analytic continuation: [Still ’89, B. ’07] Wavenumber explicit hp -estimates: [Moiola/Hiptmair/Perugia ’11]
Assembly of discrete problems Least-Squares Method Define Functional � � | [ ∇ v ] ( x ) | 2 ds + k 2 | [ v ] | 2 ds J ( v ) = Γ i ∩ Γ j i < j + least-squares for boundary conditions Discrete Solution: ˆ v := arg min v ∈V J ( v ) [Stojek ’98, Monk/Wang ’99, Barnett/B. ’10]
Scattering from a snowflake 1 0.5 0 −0.5 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 k 50 100 200 500 Time with single core Matlab code: time 8s 15s 44s 7m Least-Squares residual: ≈ 10 − 8 Exponential convergence in number of basis functions [Barnett/B. ’10]
A multiply connected domain k = 100, 30 seconds for setup and solution. Residual ≈ 2 × 10 − 5
A general plane wave DG framework [Gittelson/Hiptmair/Perugia ’09] − ∆ u − k 2 u = 0 Integrate by parts on element K with test function v : � � � ∇ u ∇ vdV − k 2 uvdV − ∇ u · nvdS = 0 ∂ K K K Further integration by parts yields: � � � ( − ∆ v − k 2 v ) udV + u · ∇ v · ndS − ∇ u · nvdS = 0 K ∂ K ∂ K • Zero for Trefftz basis functions • Replace by numerical fluxes
An example • Around 1.5 wavelengths per element • 1600 elements • 30 plane wave directions in each element • 48000 unknowns Error estimates: � u − u ( h ) � L 2 (Ω) = O ( h m + 1 ) 2 m + 1 directions (2D), ( m + 1 ) 2 directions (3D) • Framework for general Helmholtz problems • Reduction of number of unknowns compared to standard finite elements • Preconditioning difficult
Pro’s and Con’s of Trefftz methods Advantages • Highly effective up to mid-frequency regime if analytic information available [Barnett/B. ’08, Barnett/B. ’10] • Degrees of freedom recuded by constant factor compared with standard low-order FEM Disadvantages • Preconditioning difficult. Conditioning is severe problem [Huttunen/Monk/Kaipio ’02] • Restricted to homogeneous media problems • No experience whether Trefftz can beat optimized high order spectral methods • Boundary elements often better choice for waves in homogeneous media ( → BEM++ software project in development) Can we make plane wave methods fit for high-frequency problems in inhomogeneous media?
Polynomially modulated plane waves Problem: Standard plane wave methods not suitable for problems with varying speed of sound Higher order convergence only for solutions of − ∆ u − k 2 u = 0. Assume c ( x ) slowly varying on element K : c | K ≈ c K n K � p j ( x ) e i ω c K d j · x u ( x ) ≈ j = 1 Adapted DG formulation � � 1 � � ∇ u ·∇ vdx − ω 2 c ( x ) 2 uvdx + (ˆ u − u ) ∇ v · nds − ik σ · nvds = 0 ˆ K K ∂ K ∂ K
A waveguide Waveguide c ( x )
Results � u − u ( h ) � L 2 p dofs/element � u � L 2 1 30 5.7 2 60 2.3E-2 3 100 1.2E-3 ω = 10 10 directions per element
Wave travelling through a quadratic bubble • ω = 40 • 15 plane waves per element in homogeneous part • 15 modulated plane waves with degree 2 pol. in inhomogeneous part • Up to two wavelengths per element inside bubble ( 2 − r 2 ) − 1 , � r ≤ 1 c ( x ) = 1 , r > 1 Can we choose directions in a better way?
Elements from high-frequency asymptotics ∞ � High-frequency ansatz: u ( x ) = e i ωφ ( x ) a j ( x )( i ω ) − j . j = 0 φ : Solution of eikonal equation 1 |∇ φ ( x ) | = c ( x ) A : Solution of transport equation 2 ∇ φ ( x ) · ∇ A ( x ) + ∆ φ ( x ) A ( x ) = 0 Survey about high-frequency asymptotics [Engquist/Runborg ’03]
Finite wavenumber solutions Substitute u = A ( x ) e i ωφ ( x ) into Helmholtz equation � � 1 |∇ φ | 2 − 1 ω ∆ A + 2 i ∇ A · ∇ φ + iA ∆ φ − ω A = 0 c 2 Use standard polynomial approximation for A A oscillatory if • Phase approximation is bad • u has multiple phase components Amplitude FEM [Giladi/Keller ’01] Use as preconditioner [Haber/MachLachlan ’11]
Detecting wave directions • Numerical microlocal analysis [Benamou/Collino/Runborg ’04] (need to know solution) • Raytracing, fast eikonal solvers • Scattering from unit square • Each element has at most two dominant directions
Subtracting out directions Can we subtract out dominant wave directions? Assume that in element K we have n α j e ikd j · x + ˜ � u ( x ) = u j = 1 u : exact solution d j : dominant directions obtained from raytracing ˜ u : error u + k 2 ˜ Have ∆˜ u = 0, → ˜ u oscillates with same wavelength as u Approximating ˜ u not easier than approximating u !
Modulated plane waves Assume diam ( K ) = h . u ( x ) = A ( x ) e ik φ ( x ) exact solution 1 Direction d , such that ∇ φ ( x 0 ) = c ( x 0 ) d Approximate p ( x ) e i c ( x 0 ) d · x ω A ( x ) e i ωφ ( x ) ≈ � � d i ω φ ( x ) − c ( x 0 ) · x p ( x ) ≈ A ( x ) e 2 ( x − x 0 ) T ∇ 2 φ ( x 0 )( x − x 0 ) T ) + h.o.t CA ( x ) e i ω ( 1 = If d � = ∇ φ ( x 0 ) have oscillations on the scale of � h 2 2 �∇ 2 φ ( x 0 ) � , sin ∠ ( d , ∇ φ ) � ω max c ( x 0 ) Even for homogeneous problems are modulated plane waves useful!
GTD Correction Geometric optics approximation only valid for λ → 0. For finite wavelengths need GTD correction [Keller ’62]: u D ≈ c e ikr r 1 / 2 Add diffraction terms to basis n K M p j ( x ) e ikd j · x + p ℓ ( x ) H ( 1 ) � � u ( x ) ≈ 0 ( k | x − y ℓ | ) j = 1 ℓ = 1
GTD corrected basis for the square • Basis fct. with pol. of degree 3 • In each element at most 2 pw + diffraction terms • Approx 1.5 wavelengths per element • k = 40
An inhomogeneous test case 3 2 [ ( x + 1 3 ) 2 − ( y + 1 3 ) 2 ] u ( x , y ) = e i ω √ 2 Have ∆ u + ω 2 ( 3 x + 1 ) 2 + ( 3 y + 1 ) 2 � � u = 0 2 • ω = 40 • wavenumber grows from k = 20 to k = 80 across domain
Performance with plane waves 7 pw directions, L 2 error decreases quadratically [Gittelson/Hiptmair/Perugia ’09]
Convergence with exact gradient directions
Convergence with perturbed gradient directions
Recommend
More recommend
