High-order quadratures for boundary integral equations: a tutorial - PowerPoint PPT Presentation

High-order quadratures for boundary integral equations: a tutorial CBMS conference on fast direct solvers 6/23/14 Alex Barnett (Dartmouth College) Slides accompanying a partly chalk talk. Certain details, references, codes, exercises: download quadrtut.zip

Representing PDE solns: potential theory ‘charge’ (source of waves) distributed along curve Γ w/ density func. Φ Φ (x,y) ω (x,y) ρ ω Single-, double-layer potentials, x ∈ R 2 ρ n y n y � v ( x ) = Γ Φ( x , y ) σ ( y ) ds y := ( S σ )( x ) y y ∂ Φ � u ( x ) = ∂n y ( x , y ) σ ( y ) ds y := ( D σ )( x ) Γ 4 H (1) Γ Γ Φ( x , y ) := Φ( x − y ) = i 0 ( ω | x − y | ) kernel is Helmholtz fundamental soln a.k.a. free space Greens func DLP SLP

Representing PDE solns: potential theory ‘charge’ (source of waves) distributed along curve Γ w/ density func. Φ Φ (x,y) ω (x,y) ρ ω Single-, double-layer potentials, x ∈ R 2 ρ n y n y � v ( x ) = Γ Φ( x , y ) σ ( y ) ds y := ( S σ )( x ) y y ∂ Φ � u ( x ) = ∂n y ( x , y ) σ ( y ) ds y := ( D σ )( x ) Γ 4 H (1) Γ Γ Φ( x , y ) := Φ( x − y ) = i 0 ( ω | x − y | ) kernel is Helmholtz fundamental soln a.k.a. free space Greens func DLP SLP Jump relations: limit as x → Γ can depend on which side ( ± ): v ± = Sσ S, D : bdry integral ops w/ above kernels, no jump u ± = Dσ ± 1 2 σ smoothing, bounded L 2 (Γ) → H 1 (Γ) jump

Representing PDE solns: potential theory ‘charge’ (source of waves) distributed along curve Γ w/ density func. Φ Φ (x,y) ω (x,y) ρ ω Single-, double-layer potentials, x ∈ R 2 ρ n y n y � v ( x ) = Γ Φ( x , y ) σ ( y ) ds y := ( S σ )( x ) y y ∂ Φ � u ( x ) = ∂n y ( x , y ) σ ( y ) ds y := ( D σ )( x ) Γ 4 H (1) Γ Γ Φ( x , y ) := Φ( x − y ) = i 0 ( ω | x − y | ) kernel is Helmholtz fundamental soln a.k.a. free space Greens func DLP SLP Jump relations: limit as x → Γ can depend on which side ( ± ): v ± = Sσ S, D : bdry integral ops w/ above kernels, no jump u ± = Dσ ± 1 2 σ smoothing, bounded L 2 (Γ) → H 1 (Γ) jump • From now fix Γ = ∂ Ω i.e. densities live on obstacle boundary

Underlying quadrature schemes in 2D periodic composite trapezoid Gaussian panels rule err O ( N − 2 p ) if f, ∂ Ω ∈ C 2 p err O ( e − αN ) if analytic f , ∂ Ω adaptivity, corner refinement vesicles, smooth bodies

Classification of log singular schemes in 2D 4 sin 2 s − t kernel: � � K ( s, t ) = K 1 ( s, t ) log + K 2 ( s, t ) 2 split into K 1 , K 2 explicit split into K 1 , K 2 unknown global Kress ’91: prod. quadr. Kapur–Rokhlin ’97: corr. weights (PTR) Alpert ’99: aux. nodes but not FMM QBX ’12: local exp. for PDE panel-based Helsing ’08: Gen. Gauss. Kolm–Rokhlin: (Gauss–L) C contour integr. sets of aux. nodes QBX ’12 : local exp. for PDE • explicit split: more analytic info ⇒ gains efficiency • unknown split: useful for new kernels (eg axisymmetric)

Potential evaluation close to boundary 2D interior Laplace ( k = 0 ) ∂ Ω param by Z ( s ) , s ∈ [0 , 2 π ) σ ( y ) v ( z ) = i � say want eval. u = D σ u = Re v , z − ydy 2 π ∂ Ω

Potential evaluation close to boundary 2D interior Laplace ( k = 0 ) ∂ Ω param by Z ( s ) , s ∈ [0 , 2 π ) σ ( y ) v ( z ) = i � say want eval. u = D σ u = Re v , z − ydy 2 π ∂ Ω “ 5 h -rule”: target z must be 5 h from ∂ Ω to be accurate Eg use PTR. u convergence at z: log evaluation error in due to quadrature with N nodes: 10 0 0 −5 10 −5 −5 error at z −10 10 −10 −10 z z −15 −15 −15 10 N = 60 N = 120 100 200 300 N • exponential convergence, but rate arbitrarily slow as z → ∂ Ω Thm (B ’12) : rate = dist. of Z − 1 ( z ) from real axis in complex s plane

Quadrature By eXpansion (QBX) (B ’11) (Klöckner-B-Greengard-O’Neil ’12) σ , v | ∂ Ω analytic ⇒ v extends analytically some dist. outside Ω

Quadrature By eXpansion (QBX) (B ’11) (Klöckner-B-Greengard-O’Neil ’12) σ , v | ∂ Ω analytic ⇒ v extends analytically some dist. outside Ω v ( z ) = � ∞ n =0 c n ( z − z 0 ) Taylor exp. rad. of conv. ρ takes you beyond ∂ Ω u 0 log err in 10 −2 −4 Ω −6 z 0 −8 −10 −12 −14 −16 4 2 0 −2 c (n=8) integrand for n −4 −6 1 2 3 4 5 6 s

Quadrature By eXpansion (QBX) (B ’11) (Klöckner-B-Greengard-O’Neil ’12) σ , v | ∂ Ω analytic ⇒ v extends analytically some dist. outside Ω v ( z ) = � ∞ n =0 c n ( z − z 0 ) Taylor exp. rad. of conv. ρ takes you beyond ∂ Ω • pick center z 0 about 2 . 5 h from ∂ Ω u 0 log err in 10 −2 • eval. P ( ≈ 10 ) terms via Cauchy, −4 Ω −6 c n = v ( n ) ( z 0 ) = i � σ ( y ) z 0 ( z − y ) n +1 dy −8 n ! 2 π −10 ∂ Ω −12 integrand more osc. ⇒ need βN nodes, β ≈ 4 −14 interpolate σ from original N −16 4 2 • eval. Taylor exp. in | z − z 0 | ≤ R < ρ 0 −2 c (n=8) integrand for n −4 −6 • repeat for z 0 ’s all around ∂ Ω 1 2 3 4 5 6 s

Quadrature By eXpansion (QBX) (B ’11) (Klöckner-B-Greengard-O’Neil ’12) σ , v | ∂ Ω analytic ⇒ v extends analytically some dist. outside Ω v ( z ) = � ∞ n =0 c n ( z − z 0 ) Taylor exp. rad. of conv. ρ takes you beyond ∂ Ω • pick center z 0 about 2 . 5 h from ∂ Ω u 0 log err in 10 −2 • eval. P ( ≈ 10 ) terms via Cauchy, −4 Ω −6 c n = v ( n ) ( z 0 ) = i � σ ( y ) z 0 ( z − y ) n +1 dy −8 n ! 2 π −10 ∂ Ω −12 integrand more osc. ⇒ need βN nodes, β ≈ 4 −14 interpolate σ from original N −16 4 2 • eval. Taylor exp. in | z − z 0 | ≤ R < ρ 0 −2 c (n=8) integrand for n −4 −6 • repeat for z 0 ’s all around ∂ Ω 1 2 3 4 5 6 s � P + Cp � R � Cβ � P e − Cβ asymp. exponential conv. in P , β Thm. (B ’12) err ≤ C ρ P

Quadrature By eXpansion (QBX) (B ’11) (Klöckner-B-Greengard-O’Neil ’12) σ , v | ∂ Ω analytic ⇒ v extends analytically some dist. outside Ω v ( z ) = � ∞ n =0 c n ( z − z 0 ) Taylor exp. rad. of conv. ρ takes you beyond ∂ Ω • pick center z 0 about 2 . 5 h from ∂ Ω u 0 log err in 10 −2 • eval. P ( ≈ 10 ) terms via Cauchy, −4 Ω −6 c n = v ( n ) ( z 0 ) = i � σ ( y ) z 0 ( z − y ) n +1 dy −8 n ! 2 π −10 ∂ Ω −12 integrand more osc. ⇒ need βN nodes, β ≈ 4 −14 interpolate σ from original N −16 4 2 • eval. Taylor exp. in | z − z 0 | ≤ R < ρ 0 −2 c (n=8) integrand for n −4 −6 • repeat for z 0 ’s all around ∂ Ω 1 2 3 4 5 6 s | n | <P c n J n ( kr ) e inθ • Taylor → local expansion � Helmholtz ( k> 0 ): • Cauchy → Graf’s addition theorem for Bessels

QBX, 2D, high- k close eval. for Helmholtz 100 λ diameter 700 λ perimeter underlying Kress, N =9000 unknowns fill + solve 90 sec QBX eval in 30 sec ( 2 × 10 5 pts) rel. error < 10 − 11 (B ’12)

Local vs global QBX; same scheme in 3D Local: use QBX to fill self and near panel matrix blocks, sparse O ( N ) – far via plain rule; err O ( h p + ǫ ) where ǫ fixed, controlled by p, P, β . ie not formally convergent; needs P high to push to ǫ = O ( ǫ mach ) Global: use QBX with all of ∂ Ω contrib to expansion at each center – kills the ǫ , allows lower P (for engs. apps.), do all via FMM (FDS?)

Local vs global QBX; same scheme in 3D Local: use QBX to fill self and near panel matrix blocks, sparse O ( N ) – far via plain rule; err O ( h p + ǫ ) where ǫ fixed, controlled by p, P, β . ie not formally convergent; needs P high to push to ǫ = O ( ǫ mach ) Global: use QBX with all of ∂ Ω contrib to expansion at each center – kills the ǫ , allows lower P (for engs. apps.), do all via FMM (FDS?) fine source mesh for QBX centers for this panel self and neighboring panels 3D: panels p × p Gauss nodes panel of targets (p p, eg p=8) u ( r, θ, φ ) = Local expansion n � � c nm j n ( kr ) Y m n ( θ, φ ) m = − n | n |≤ P spherical harmonic addn thm • ( P +1) th order proven for σ ∈ W P +3+ ǫ, 2 (Epstein–Greengard–Klöckner ’12)

QBX: 3D high-freq. torus scattering result Dirichlet BC (sound-soft acous- tics) 30 λ diameter N ≈ 145000 q =8 , p =10 , β =4 . 5 QBX quad 1.2 hr GMRES+FMM 1 hr laptop (4-core i7) relative error 10 − 5 • QBX in 3D still in primitive state (Barnett–Gimbutas–Greengard, in prep.) • note FEM/FDTD at this high accuracy & freq. essentially prohibitive

QBX: 3D periodic scattering (prelim) Doubly-periodic grating of sound-soft scatterers Dirichlet obstacles d = 2 . 4 λ N = 25200 (one obstacle) p = 6 . QBX 4 min, laptop p =6 , P =8 , β =4 30 its 5 min error 10 − 5 • New periodizing scheme (Barnett–Gimbutas–Greengard, in prep.)

3D, Bremer–Gimbutas ’12: triangle auxiliary nodes Lots of precomputed nodes for various aspect triangles, kernels: local correction (self & neighbors) product grids in two parameters polar coords removes 1 /r singularity Low-frequency Helmholtz Neumann BVP: