Stable Outgoing Wave Filters for Anisotropic Waves Chris Stucchio - PowerPoint PPT Presentation

Stable Outgoing Wave Filters for Anisotropic Waves Chris Stucchio Courant Institute of Mathematical Sciences Joint work with Avy Soffer. Supported by NSF and Bevier Fellowship. Monday, July 7, 2008 1

Linear Waves • Linear wave equation: u t ( x, t ) = u ( , t ) � H � − H † ( i ∇ ) H ( i ∇ ) = Monday, July 7, 2008 2

Linear Waves • Linear wave equation: u t ( x, t ) = u ( , t ) � H � − H † ( i ∇ ) H ( i ∇ ) = • Schrodinger equation: = i ∆ H u ( x, t ) = ψ ( x, t ) Monday, July 7, 2008 3

Linear Waves • Linear wave equation: u t ( x, t ) = u ( , t ) � H � − H † ( i ∇ ) H ( i ∇ ) = • Maxwell’s equation � 0 − µ − 1 / 2 ∇ × ǫ − 1 / 2 � = H ǫ − 1 / 2 ∇ × µ − 1 / 2 0 H, √ ǫ � ( √ µ � u ( x, t ) = E ) � Monday, July 7, 2008 4

Linear Waves • Linear wave equation: u t ( x, t ) = u ( , t ) � H � − H † ( i ∇ ) H ( i ∇ ) = • Linearized Euler equation:   M ∂ x 1 − ∂ x 1 − ∂ x 2 = 0 − ∂ x 1 M ∂ x 1 H   0 − ∂ x 2 M ∂ x 1 ( x, y ) = ( p ( x, t ) , v x ( x, t ) , v y ( x, t )) Monday, July 7, 2008 5

Linear Waves • Linear wave equation: u t ( x, t ) = u ( , t ) � H � − H † ( i ∇ ) H ( i ∇ ) = • Relativistic Schrodinger Equation � = − ∆ + m 2 − m H u ( x, t ) = ψ ( x, t ) Monday, July 7, 2008 6

Linear Waves • Linear wave equation: u t ( x, t ) = u ( , t ) � H � − H † ( i ∇ ) H ( i ∇ ) = • Linear part of Benjamin-Ono equation: = | ∂ x | ∂ x H u ( x, t ) = h ( x, t ) Monday, July 7, 2008 7

Numerical Solution • Finite Differences • Finite Elements • Spectral methods I’ll stay agnostic FFT spectral methods rock. Monday, July 7, 2008 8 Mention Fundamental Complexity

Numerical Solution • Sample spacing: δ x ≤ O (2 π /k max ) [ − L, L ] N • Fundamental complexity of timestepping on O (( Lk max ) N ) Memory = O (( T max / δ t )( Lk max ) N ) Complexity = • Solution on requires careful choice of boundary conditions. R N Monday, July 7, 2008 9

Outgoing Waves are a Problem Monday, July 7, 2008 10

The Problem 1D Schrodinger Equation Monday, July 7, 2008 11

The Problem Anisotropic Maxwell Incorrect Boundaries Monday, July 7, 2008 12 Quality is a MAC issue, not an accuracy one. Looks better on Linux.

Possible Solution: Exact NRBC • Dirichlet-to-Neumann boundaries: impose exact non-reflecting boundary conditions, constructed from Green’s function to free wave. • Nonlocal in time, nonlocal on boundary • Internal solver restricted (no Fourier spectral methods) • Geometry restricted • Majda-Engquist, Bayliss-Turkell, Hagstrom, Greengard, Grote, ... Monday, July 7, 2008 13

Possible Solution: Perfectly Matched Layers • Extend with absorbing layer PML • Dissipation inside layer • Must be Perfectly Matched to Region avoid reflection at the interface. of • Equivalent to complex scaling Interest Monday, July 7, 2008 14

Possible Solution: Perfectly Matched Layers • Complex scaling for Wave equation: e zA He − zA H �→ A x · i ∇ + i ∇ · x = • PML (Conjugate Operator) for general linear waves: e zA He − zA H �→ = x · v g ( i ∇ ) + v g ( i ∇ ) · x A Monday, July 7, 2008 15

Complex scaling is easy = x · i ∇ + i ∇ · x A e zA = Dilation(z) • Change coordinates • Make layer perfectly matched Picture from Becache, Fauqueux, k 1 v g, 1 ( k ) ≥ 0 • Stable if Joly, JCP 188 (2003) 399–433. Monday, July 7, 2008 16

PML Instability • PML unstable for some anisotropic waves (Becache, Fauqueux, Joly, 2003). Pictures from Becache, Fauqueux, Joly, JCP 188 (2003) 399–433. Monday, July 7, 2008 17 Picture and graph are for di fg erent equations

Conjugate operators are hard = x · v g ( i ∇ ) + v g ( i ∇ ) · x A e zA = ? Monday, July 7, 2008 18 If stability condition not satisfied, complex scaling shifts spectrum the wrong way.

Phase Space Filters • Identify outgoing waves • Filter them off • Nothing hits the boundary Monday, July 7, 2008 19

Outgoing waves Monday, July 7, 2008 20

Outgoing Waves, Schrodinger Equation • 1D Schrodinger Equation � − x 2 e ivx � ψ 0 ( x ) = √ σ exp 2 σ 2 � − ( x − vt ) 2 e ivx � ψ ( x, t ) = exp � 2 σ 2 (1 + it/ σ ) σ + it/ σ x = vt width = σ + t/ σ • Center of mass at , Monday, July 7, 2008 21

Outgoing Waves, Schrodinger Equation • Outgoing wave e + ivx e − ( x − L ) 2 / σ 2 ψ 0 ( x ) = Trajectory = L + vt • Incoming wave e − ivx e − ( x − L ) 2 / σ 2 ψ 0 ( x ) = Trajectory = L − vt Monday, July 7, 2008 22

Outgoing Waves, Schrodinger Equation • Outgoing wave e + ivx e − ( x − L ) 2 / σ 2 ψ 0 ( x ) = Trajectory = L + vt • Incoming wave e − ivx e − ( x − L ) 2 / σ 2 ψ 0 ( x ) = Trajectory = L − vt Monday, July 7, 2008 23

Outgoing Waves, Schrodinger Equation • Outgoing wave e + ivx e − ( x − L ) 2 / σ 2 ψ 0 ( x ) = Trajectory = L + vt • Incoming wave e − ivx e − ( x − L ) 2 / σ 2 ψ 0 ( x ) = Trajectory = L − vt Monday, July 7, 2008 24

Outgoing Waves, Schrodinger Equation • Mixed wave: e − ivx e − ( x − L ) 2 / σ 2 + e + ivx e − ( x − L ) 2 / σ 2 ψ 0 ( x ) = Monday, July 7, 2008 25

Outgoing Waves, Schrodinger Equation • Mixed wave: e − ivx e − ( x − L ) 2 / σ 2 ψ 0 ( x ) = Incoming wave Outgoing wave Monday, July 7, 2008 26

Outgoing Waves, Schrodinger Equation • Mixed wave: e − ivx e − ( x − L ) 2 / σ 2 ψ 0 ( x ) = Incoming wave Outgoing wave Monday, July 7, 2008 26

Outgoing Waves, Schrodinger Equation • Mixed wave: e − ivx e − ( x − L ) 2 / σ 2 ψ 0 ( x ) = Incoming wave Problem solved! Monday, July 7, 2008 27

It really is that easy • Windowed Fourier Transform: � � ψ a,b e ibk 0 x g ( x − ax 0 ) ψ ( x ) = a ∈ Z b ∈ Z g ( x ) = e − x 2 / σ 2 • Outgoing waves: ax 0 > L The Wavelet Transform, Time Frequency σ − 1 Localization and Signal Analysis , Ingrid bk 0 > Daubechies, IEEE Trans. Info. Theory, Vol 36 5 1990 Monday, July 7, 2008 28

Quantum Phase Space • Quantum phase space is set of points , x a ( x, k ) ∈ R N × R N position and k a frequency. • Heisenberg Uncertainty principle: localizing on region of volume causes O (2 π ln( ǫ )) error . ǫ • A function is localized near a point if it is localized in ( x 0 , k 0 ) position near and it’s x 0 Fourier transform is localized near . k 0 Monday, July 7, 2008 29

Outgoing Waves, Schrodinger Equation • Ambiguous waves ψ 0 ( x ) = e i 0 x e − ( x − L ) 2 / σ 2 • Spreads in both directions Issue can be resolved. Monday, July 7, 2008 30

Phase space filters Monday, July 7, 2008 31

Outgoing ax 0 > L Phase Space Filters Waves: σ − 1 bk 0 > Monday, July 7, 2008 32

Outgoing ax 0 > L Phase Space Filters Waves: σ − 1 bk 0 > Monday, July 7, 2008 32

Outgoing ax 0 > L Phase Space Filters Waves: σ − 1 bk 0 > Monday, July 7, 2008 32

A simpler version • We want rightward moving waves near the boundary . • Extend computational domain [ − L − w, L + w ] N • Localize in boundary layer 1( x > L )1( k > k min )1( x > L ) = O + Monday, July 7, 2008 33 Step function built out of erfc(x) functions, as smooth as possible without overflowing.

How does it work? • Take wave comprised of incoming and outgoing waves, plus interior waves. e ivx g ( x − L − 1) + e − ivx g ( x − L − 1) ψ 0 ( x ) ≈ + interior waves Monday, July 7, 2008 34

How does it work? • Take wave comprised of incoming and outgoing waves, plus interior waves. e ivx g ( x − L − 1) + e − ivx g ( x − L − 1) ψ 0 ( x ) ≈ + interior waves Our target Monday, July 7, 2008 35

How does it work? • Take wave comprised of incoming and outgoing waves, plus interior waves. e ivx g ( x − L − 1) + e − ivx g ( x − L − 1) ψ 0 ( x ) ≈ + interior waves • We don’t care about interior waves e ivx g ( x − L − 1) + e − ivx g ( x − L − 1) 1( x > L ) ψ 0 ( x ) ≈ + 0 Monday, July 7, 2008 36

How does it work? • Take wave comprised of incoming and outgoing waves, plus interior waves. e ivx g ( x − L − 1) + e − ivx g ( x − L − 1) ψ 0 ( x ) ≈ + interior waves • We don’t care about interior waves e ivx g ( x − L − 1) + e − ivx g ( x − L − 1) 1( x > L ) ψ 0 ( x ) ≈ + 0 • Or incoming waves e ivx g ( x − L − 1) + 0 1( k > k min )1( x > L ) ψ 0 ( x ) ≈ Monday, July 7, 2008 37

How does it work? • Or incoming waves e ivx g ( x − L − 1) + 0 1( k > k min )1( x > L ) ψ 0 ( x ) ≈ • Symmetry is always good (for stability, etc): e ivx g ( x − L − 1) 1( x > L )1( k > k min )1( x > L ) ψ 0 ( x ) ≈ Monday, July 7, 2008 38