Hybrid (DG) Methods for the Helmholtz Equation Joachim Sch oberl - PowerPoint PPT Presentation

Hybrid (DG) Methods for the Helmholtz Equation Joachim Sch¨ oberl Computational Mathematics in Engineering Institute for Analysis and Scientific Computing Vienna University of Technology Contributions by A. Sinwel (Pechstein), P. Monk M. Huber, A. Hannukainen, G. Kitzler RICAM Workshop, Nov 2011, Linz Joachim Sch¨ oberl Page 1

Helmholtz Equation Heterogeneous Helmholtz Equation ∆ u + n ( x ) 2 ω 2 u = f Impedance trace boundary condition ∂u ∂n + iωu = iωg in for convenience. better: PML, pole-condition, ... Weak form: Find u ∈ H 1 s.t. � � � ∇ u ∇ v − n 2 ω 2 uv + iω ∀ v ∈ H 1 uv = gv Ω ∂ Ω ∂ Ω Joachim Sch¨ oberl Page 2

Weak and dual formulation Mixed formulation. Introduce velocity σ := 1 iω ∇ u : ∇ u − iωσ = 0 1 div σ − iωn 2 u = iωf σ n + nu = g in Dual formulation. Eliminate u : ∇ 1 n 2 div σ + ω 2 σ = ∇ 1 n 2 f Weak form: � � � � � i 1 1 n 2 div σ div τ − iω στ + nσ n τ n = f div τ + g in τ n ω ∂ Ω ∂ Ω Joachim Sch¨ oberl Page 3

Function space H (div) Dual formulation requires the function space H (div) = { τ ∈ [ L 2 ] d : div τ ∈ L 2 } has well-defined normal-trace operator Raviart-Thomas Finite Elements: a ∈ [ P k ] d , b ∈ P k } V RT k = { � a + b� x : � There holds [ P k ] d ⊂ V RT k ⊂ [ P k +1 ] d , div V RT k = P k , tr n,F V RT k ∈ P k ( F ) Degrees of freedom are • moments of normal traces on facets F of order k • moments on the element T of order k − 1 dofs ensure continuity of normal-trace Joachim Sch¨ oberl Page 4

Hybridization of RT -elements Arnold-Brezzi hybridization technique: Do not incorporate normal-continuity σ | T 1 · n 1 = − σ | T 2 · n 2 into the finite element space, but enforce it by additional equations: � [ σ n ] v F = 0 ∀ v ∈ P k ( F ) F RT k , u F ∈ P k ( F ) such that Hybrid system: Find σ ∈ V dc � � � � � i ∂T τ n u F T div σ div τ − iω T στ + = 0 ∀ τ T T ω � � � � ∂T σ n v F ∂ Ω u F v V ∂ Ω g in v F ∀ v F + = T u F is the primal variable on the facet. Wish to eliminate σ from first equation (small local problems !), BUT local Dirichlet - problems are unstable for hω � 1 Joachim Sch¨ oberl Page 5

Hybridization of Helmholtz equation F u σ F Add a second Lagrange multiplier [Monk+Sinwel+JS 11] n σ F σ n F := σ | F · n F σ � ∂T ( σ n − σ F n )( τ n − τ F and stabilize by adding 0 = n ) � � � � � � ∂T τ n u F − τ n σ F i T div σ div τ − iω T στ + ∂T σ n τ n + = 0 ∀ τ n T T ω � ∂T σ n v F − σ n τ F � � � � � ∂T σ F n τ F ∂ Ω u F v F ∂ Ω g in v F ∀ v F + n + = n T T The element-variable σ can now be eliminated from the stable local equation (Robin - b.c.): � � � � i ( σ F n − u F ) τ n div σ div v − iω στ + σ n τ n = ω T ∂T ∂T Joachim Sch¨ oberl Page 6

Analysis • Unique discrete solution and convergence rate estimates by usual duality technique ( h sufficiently small) • Qualitative property: Impedance trace isometry • Local solvability (by Melenk test-functions τ = x div σ ) Joachim Sch¨ oberl Page 7

Impedance trace isometry Element equation reads as � � � � i ( σ F n − u F div σ div v − iω στ + σ n τ n = ) τ n ω � �� � T ∂T ∂T g in define element-wise out-going impedance trace g out ∈ P k ( F ) via � � � � i div σ div v − iω στ − σ n τ n = g out τ n ω T ∂T ∂T � � � Choose τ = σ and take real parts: � σ n � 2 ∂T = R g in σ n So we get � ∂T = � g in � 2 − 4 R g in σ n + 4 � g in � 2 = � g in � 2 � g out � 2 ∂T = � g in − 2 σ n � 2 ∂T The facet-equations can be rephrased as (for neighbour elements ˜ T ). ˜ T ) g in = g out ( σ Joachim Sch¨ oberl Page 8

Convergence plots Comparing Hybrid-Helmholtz solution with conforming H 1 -FEM. n , u F ∈ P k ( F ) . Postprocessing to ˜ σ ∈ RT k , σ F u ∈ P k +1 ( T ) . L-2 error u 10 RT3 - postproc std FEM 4 1 0.1 0.01 0.001 err 0.0001 1e-05 1e-06 1e-07 1e-08 0.01 0.1 1 h Hybrid-Helmholtz has about twice the global dofs of standard FEM. Joachim Sch¨ oberl Page 9

How to solve ???? Want to iterate for facet-variables only. Idea: In the wave-regime the wave-relaxation ˆ T ) g in := g out ( σ converges well. But, does not work for the elliptic regime when hω is small. Try to combine with an elliptic preconditioner. Joachim Sch¨ oberl Page 10

Facet dofs Facet - dofs are Dirichlet and Neumann data u F , σ F n , or left- and right-going traces g l , g r transformation � � � � � � g l u F 1 1 = g r σ F 1 − 1 n holds point-wise since all variables are in P k ( F ) Joachim Sch¨ oberl Page 11

Motivation: 1D system The 1D system (with exact local problems) reads as � � g l � � g l � � g l � � � � � � 0 γ 0 1 0 0 j − 1 j j +1 + + g r g r g r 0 0 1 0 γ 0 j +1 , j − 1 j left and right-going plane waves with phase-shift γ . left (right) going Gauss-Seidel for g l ( g r ) is an exact solve. Matrix is not diagonal dominant. But, exchanging rows pairwise gives a diagonal-dominant, non-symmetric matrix. Krylov-space methods with most diagonal-like preconditioners will converge. Joachim Sch¨ oberl Page 12

Balancing domain decomposition methods Domain decomposition, Ω = ∪ Ω i , distribute elements Sub-assemble fe matrices with sub-domain matrix A i � A = A i Balancing DD Preconditioner: � A − 1 = ˆ R i A − 1 i R T i Elliptic case: local Neumann problems are singular ! Solution: BDD with constraints ⇒ BDDC (Dohrmann, Widlund, ...)   − 1 A cc A c 1 · · · A cm A 1 c A 11   R T R   . ... . .   A mc A mm Analysis for the elliptic case with log α condition numbers. Joachim Sch¨ oberl Page 13

Heuristics: BDDC for the Hybrid Helmholtz system Bilinearform � � � i 1 � σ n v F + τ n u F + B ( σ, σ F , u F ; τ, τ F , v F ) = div σ div τ − iω στ + ω ∂T T � � n v F + τ F ( σ n − σ F n )( τ n − τ F σ F n u F + n ) − γ ∂T ∂T The new blue term is consistent and is motivated by the g l − g r row-exchange. Apply BDDC preconditioner with mean-value facet-dof constraints Good values for γ are found by experiments. Joachim Sch¨ oberl Page 14

Iteration numbers with coarse grid, 4 domains with coarse grid, 16 domains L/λ 2 8 32 L/λ 2 8 32 p=2 21 25 - p=2 23 32 - p=4 27 30 42 p=4 33 41 76 p=8 38 46 41 p=8 42 57 71 p=16 59 73 60 p=16 65 87 100 no coarse grid, 4 domains no coarse grid, 16 domains L/λ 2 8 32 L/λ 2 8 32 p=2 51 40 - p=2 85 80 - p=4 64 48 40 p=4 109 88 101 p=8 96 61 43 p=8 148 110 89 p=16 132 84 63 p=16 214 159 123 Joachim Sch¨ oberl Page 15

Infrastructure • general purpose hp - Finite Element Code Netgen/NGSolve C++, open source MPI - parallel • Dell R 910 Server 4 Xeon E7 CPUs 10 core @ 2.2 GHz 512 GB RAM Joachim Sch¨ oberl Page 16

Diffraction from a grating Sphere with D = 40 λ , 127 k elements, h ≈ λ , p = 5 , 39 M dofs (corresponds to 9 . 4 M primal dofs) 78 sub-domains / processes, no coarse grid T ass = 9 m , T pre = 12 m , T solve = 21 m , 156 its Joachim Sch¨ oberl Page 17

Computational Optics factor out wavy part lense, n = 2 u = e ik · x ˜ u � � ω simple guess: k = n 0 better: ray-tracing, Eukonal equation wavy u smooth ˜ u Joachim Sch¨ oberl Page 18