Modulated plane wave methods for Helmholtz problems in heterogeneous - - PowerPoint PPT Presentation

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

∆u + ω c(x) 2 u = k = ω c

• Trefftz methods for

homogeneous media

• Polynomially modulated plane

waves in inhomogeneous media

• 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 := −∆ − k2). Ω :=

• i

Ki: Decomposition of domain Ω into elements Ωi. Local approximation spaces Vi := {

n

• j=1

α(i)

j φj : Lφj = 0 in Ki}.

• Local vs. global methods
• Choices of basis functions
• Assembly of discrete problem

Trefftz basis functions

Plane Waves φj = eikdj·x, dj = 1 dj = cos θj sin θj

• , θj = 2π(j − 1)

N Fourier-Bessel functions φj(r, θ) = Jj(kr)eijθ j = −n, . . . , n

Trefftz basis functions. . .

Fundamental Solutions φj(x) = H(1)

0 (k|x − yj|)

yj: 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

Vekua operator (2d version): Vj[φ](x) = φ(x) + 1 Mj(x, t)φ(tx)dt M1(x, t) = − k|x| 2 √ 1 − t J1(k|x| √ 1 − t) M2(x, t) = −i k|x| 2 √ t √ 1 − t J1(ik|x|

• t(1 − t))

Properties:

1 V1[V2[φ]] = V2[V1[φ]] = φ. 2 If φ is harmonic, then

∆V1[φ] + k2V1[φ] = 0.

3 If ∆u + k2u = 0, then

∆V2[u] = 0 [Eisenstat ’74, Still ’89, Melenk ’99, B. ’07, Moiola/Hiptmair/Perugia ’11]

Approximation properties. . .

Generalized harmonic polynomials P(r, θ) =

n

• j=−n

ajr |j|eijφ V1[P](r, θ) =

n

• j=−n

aj|j|! 2 k |l| eijθJ|j|(kr) u −

n

• j=1

αjeidj·x ≤ V1φ − P + V1[P] −

n

• j=1

αjeidj·x

• 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 J(v) =

• i<j
• Γi∩Γj

| [∇v] (x)|2ds + k2| [v] |2ds + 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

−0.2 0.2 0.4 0.6 0.8 1 1.2 −0.5 0.5 1

Time with single core Matlab code: k 50 100 200 500 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 − k2u = 0 Integrate by parts on element K with test function v:

• K

∇u∇vdV − k2

• K

uvdV −

• ∂K

∇u · nvdS = 0 Further integration by parts yields:

• K

(−∆v − k2v)udV +

• ∂K

u · ∇v · ndS −

• ∂K

∇u · nvdS = 0

• 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)L2(Ω) = O(hm+1) 2m + 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 − k2u = 0. Assume c(x) slowly varying on element K: c|K ≈ cK u(x) ≈

nK

• j=1

pj(x)eiωcK dj·x Adapted DG formulation

• K

∇u·∇vdx−ω2

• K

1 c(x)2 uvdx+

• ∂K

(ˆ u−u)∇v · nds−ik

• ∂K

ˆ σ·nvds = 0

A waveguide

Waveguide c(x)

Results

p dofs/element

u−u(h)L2 uL2

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

c(x) =

• (2 − r 2)−1,

r ≤ 1 1, r > 1 Can we choose directions in a better way?

Elements from high-frequency asymptotics

High-frequency ansatz: u(x) = eiωφ(x)

∞

• j=0

aj(x)(iω)−j. φ: Solution of eikonal equation |∇φ(x)| = 1 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)eiωφ(x) into Helmholtz equation 1 ω ∆A + 2i∇A · ∇φ + iA∆φ − ω

• |∇φ|2 − 1

c2

• A = 0

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 u(x) =

n

• j=1

αjeikdj·x + ˜ u u: exact solution dj: dominant directions obtained from raytracing ˜ u: error Have ∆˜ u + k2˜ 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)eikφ(x) exact solution Direction d, such that ∇φ(x0) = 1 c(x0)d Approximate p(x)ei

ω c(x0) d·x

≈ A(x)eiωφ(x) p(x) ≈ A(x)e

iω

• φ(x)−

d c(x0) ·x

• =

CA(x)eiω( 1

2 (x−x0)T ∇2φ(x0)(x−x0)T) + h.o.t

If d = ∇φ(x0) have oscillations on the scale of ω max h2 2 ∇2φ(x0), sin ∠(d, ∇φ) c(x0)

• 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]: uD ≈ c eikr r 1/2 Add diffraction terms to basis u(x) ≈

nK

• j=1

pj(x)eikdj·x +

M

• ℓ=1

pℓ(x)H(1)

0 (k|x − yℓ|)

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

u(x, y) = eiω

3 2 √ 2[(x+ 1 3 )2−(y+ 1 3 )2]

Have ∆u + ω2 2

• (3x + 1)2 + (3y + 1)2

u = 0

• ω = 40
• wavenumber grows from

k = 20 to k = 80 across domain

Performance with plane waves

7 pw directions, L2 error decreases quadratically [Gittelson/Hiptmair/Perugia ’09]

Convergence with exact gradient directions

Convergence with perturbed gradient directions

Raytracing

Define Hamiltonian H(x, p) = c(x)|p| Hamiltonian constant along bicharacteristics defined by dx dt = ∇pH(x, p) = c(x) p |p| dp dt = −∇xH(x, p) = −|p|∇c(x) Look for solutions of H(x, p) = 1. x(t) : Ray p(t) : Slowness vector For homogeneous material rays are simple straight lines

A bubble example

Slowness function from [Engquist/Ying ’11]

Wavefront tracking [Vinje/Iversen/Gjøystdal ’92]

Summary

• Trefftz methods show good performance for homogeneous

problems

• Can we beat well optimized hp-finite elements or boundary

elements?

• For inhomogeneous problems can use modulated plane waves
• Competitive against standard methods only if high-frequency

information can be use Goal: Analyze Helmholtz problems using high-frequency asymptotics. Then refine solution with high-frequency finite elements.