Exponentially Convergent Sparse Discretizations and Application to - - PowerPoint PPT Presentation

Exponentially Convergent Sparse Discretizations

and Application to Near Surface Geophysics

North Carolina State University

November 9, 2017

Murthy N. Guddati

Outline

 Part 1: Impedance Preserving Discretization  Part 2: Absorbing Boundary Conditions (1-sided DtN map)

 Joint work with Tassoulas, Druksin, Lim, Zahid, Savadatti, Thirunavukkarasu

 Part 3: Complex-length FEM for finite domains (2-sided DtN map)

 Joint work with Druskin, Vaziri Astaneh

 Part 4: Inversion for Near-surface Geophysics

 Joint work with Vaziri Astaneh

2

• Part1. Impedance Preserving Discretization

Model Problem

 3D wave equation in free space  Fourier transform in t, y, z, with  Exact solution:  . is a plane/evanescent wave

2 2 2 2 2 2 2 2 2

1 u u u u x y z c t             

y z

ik y ik z i t

u Ue

  



2 2 2

U k U x     

2 2 2 2

is complex valued

y z

k k k c            , where is the horizontal wavenumber

ikx ikx

U Ae Be k



 

( )

y z

i kx k y k z t ikx

U e u e

   

  

4

Finite Element Solution on a Uniform Grid in x

 FE discretization of:  Element contribution matrix with uniform element size of h:  Assembly results in the difference equation:

2 2 2

U k U x     

2

1 1 1 1 1 3 6 1 1 1 1 6 3

elem

k h A A h B B                                k

2 2 2 2

1 1 3 1 1 6 k h A h k h B h                 

1 1

2

j j j

BU AU BU

 

  

5

Changing Mesh Size: Reflections

 A simple analysis using two uniform meshes with different element sizes (h, H), but the same material  What happens when a right propagating wave hits the interface?  Exact solution – just passes through  Finite element solution – reflections due to impedance mismatch

: discrete impedance of left domain : discrete impedance of right domain

h H h H H h

Z Z Z R Z Z Z   

6

Computing Discrete Impedance (Half-space Stiffness)

 Basic idea: discrete half-space + finite element = discrete half-space  . depends on element size, resulting impedance mismatch when the element size changes, resulting in reflections

2 2 2 1 1 h h h h h

A B A Z B U U Z U A Z B B A Z B A Z U U                                             

 

2 2 2

2 1 1

h

kh B ik Z A     

2 2 2 2

1 1 3 1 1 6 k h A h k h B h                 

h

Z

Error term

7

Optimal Integration for Minimizing Reflection Error

 Minimize the error in impedance by using generalized integration rules  Minimize the error term by choosing  The error in impedance is completely eliminated! No more reflections  Formally valid for more general 2nd order equations (anisotropic, visco- elasticity etc., electromagnetics etc. – G, 2006, CMAME)

 

2 2 2 2

1 4

h

kh Z A B ik      

2 2 2 2 2 2

1 1 1 4 1 1 1 4 A k h h B k h h                                 

Linear elements + midpoint integration = Impedance Preserving Discretization Error term

 

 

   

8

Part 2. Absorbing Boundary Conditions Perfectly Matched Discrete Layers

Perfectly Matched Discrete Layers

…Impedance Preserving Discretization of PML

 Perfectly Matched Layers (PML) (Berenger, 1994; Chew et.al. 1995)

 Step 1: Bend the domain into complex space  Step 2: discretize PMDL domain (in complex space)  Impedance is no longer preserved; perfect matching is destroyed  Requires a large number of carefully chosen PML layers

 Impedance preserving discretization comes to the rescue!

 Impedance is preserved/matched, irrespective of element length, small, large, real, complex – Perfectly Matching Discrete Layers (PMDL)  Discretize with 3-5 complex-length linear finite elements  No discretization error, but truncation causes

• reflections. The reflection coefficient is derived as

PML Region (imaginary or complex x) Interior (real x)

Reduced reflection into the interior

2 1

1 / 2 1 / 2

j nlayer j PMDL j j

ikL R ikL

 

          



10

2ikL PML

R e 

PMDL vs PML: Effectiveness of Midpoint Integration

PMDL with 3 layers PML with 3 layers

11

 Impedance preservation property is valid for any equation that is linear and second order in space (G, CMAME, 2006)  Elastic and other complicated wave equations (G, Lim & Zahid, 2007)  Evanescent waves can be treated effectively  Padded PMDL – contains large real lengths with midpoint integration (Zahid & G, CMAME, 2006)

PMDL with 5 layers PML with 5 layers

PMDL: Some More Old Results

12

Salient Features of PMDL

 Exponential convergence  Near optimal discretization

 Optimal: need staggered grids (with Druskin et al., 2003)

 Links PML to rational ABCs

 Lindman, Engquist-Majda, Higdon and variants (e.g. CRBC)  We started this from E-M/Higdon ABCs (G, Tassoulas, 2000)  Extensions to corners is straightforward

 Additional advantage: Provides solutions to some difficult cases

 Backpropagating waves: anisotropy  PML for discrete/periodic media

13 2 1 j nlayer j j j

k k R k k

 

          



PML Region Interior

PMDL for Backpropagating Waves

Opposing signs of phase and group velocities

 Backpropagating waves grow in the PML region

 PML cannot work! (Bécache, Fauqueux and Joly, 2003)

PML result: radiation in anisotropic media Result from PMDL after the fix

Savadatti & G (2012), J Comp. Phys.

 A counter-intuitive idea: make the reflections in PML region decay faster than the growth of the incident wave

Reduced Reflection into the interior

 Works only with PMDL: needs impedance preserving discretization!

14

Anisotropic elasticity – Tilted Elliptic Case

Stable parameters Arbitrary parameters

15

Savadatti & G (2012), J Comp. Phys.

Ideal Slowness

Anisotropic elasticity – Non-elliptic Case

16

Two different coordinated materials Traditional mesh

Savadatti & G (2012), J Comp. Phys.

PMDL for Periodic Media (after discretization)

 Periodic media has internal reflections and transmissions

 Constructive interference leads to long-range propagation

 PML’s complex stretching spoils this balance and internal reflections and transmissions get mixed up!  Basic Ideas (Discrete/Periodic PMDL):

 Periodic media = Discrete vector wave equation (vector size = ndof in a cell)  Discrete vector equation = impedance preserving discretization of more complicated wave equation  Apply PMDL on the complicated wave equation results in impedance matching for periodic media

 Open problem: stability for complex problems

PML for Lattice Waves: 7% reflections w/20 PML layers Discrete PMDL: less than 1% error w/ 4 PMDL layers 17

G & Thirunavukkarasu, JCP (2009), Waves 2011

Part 3. Two-Sided DtN Map Complex-length Finite Element Method

Facilitating the Approximation of 2-Sided DtN Map

 Consider the equation:  Exact 2-sided DtN map:  By definition, exact DtN Map is impedance preserving:  Consider impedance preserving discretization of the interval:  Error in A and B would be similar since:  Approximating two-sided map reduces to approximating one-sided map  Better derivation based on Crank-Nicolson discretization of the propagator

19

2 2 2

0, (0, ) u u z L z       

=

exact

A B K B A      

2 2 2 exact

A B Z  

2 2 2

= ,

exact exact

A B K A B Z B A        

2 2 2 2 2 exact

A B Z A B    

1D Helmholtz Equation

20

z

2 2 2

u u z      

exp(-iL) exp(+iL)

z  z L 

/ / u i i u v i i v i z                          

/ v u z   

1st Order Form

Eigenvalues i   Downwardwaves:

1

L

u i u v i v z                       

21

Propagator Approximation

. . .

1

L

2

L

n

L

1 1 1 1 1

1 1 2 u u u u i v v v v i L                         

1st Order Form Crank- Nicolson

u

1

u u

L

u

Crank- Nicolson

Valid for both Downward AND upward waves!  

/ , / v u z v v i    

approximate n

P

Padé Approximant

22

Imaginary Real

. . .

1

L

2

L

n

L

 

Padé

exp( )

j j j j

d L d P d d

 

  

 



( 0, ,2 ) j n 

(2 )! ( ) 2 / ! ( )!

n j j j j

n j x L L x j n j



     



Complex-Length Finite Element Method

Complex-Length FEM: Exponential Convergence

Laplace Equation Helmholtz Equation (can’t beat Nyquist limit)

1 L  1 f 

Complex-Length FEM: Some Observations

 Exponentially convergent  Piecewise linear interpolation – sparse computation  Edges do not move (∑Lj=L) – can be combined with other types of meshes for other subdomains  Mesh is not bent outside (Re(Lj)>0)  Order of elements do not matter! – more on this later.  With refinement, and proper ordering, mesh converges to a smooth curve on the complex plane

(2 )! ( ) 2 / ! ( )!

n j j j j

n j x L L x j n j



     



 

Padé

exp( )

j j j j

d L d P d d

 

  

 



Energy Conservation and Eigenvalue Problems

 Do complex lengths lead to energy absorption, like PML?  No, due to conjugate pair of lengths – decay grows back!  2-sided DtN Map is Hermitian

2

( )    K M

2 2 2

u u z      

 Do complex lengths lead to complex eigenvalues of K with respect to M? 

• No. Eigenvalues are real

and positive!  Eigenvectors are complex (K and M are complex symmetric)

26

Element Ordering

exp( ) u c i x   

Part 4. Near Surface Geophysical Site Characterization… …using Guided Wave Inversion

Guided Wave Dispersion

Low Frequency High Frequency

Dispersion Curve

28

Spectral Analysis

FFT Experimental Dispersion Curve Phase Velocity (m/s) = Frequency (1/s) Wavenumber (1/m)

29

30

Medium Characterization

Experimental Dispersion Curve

Inverse Identification

 

2 experimantal predicted 1 N i i i

E c c



 



Iteratively Minimize

 Optimization Scheme

 Gradient Based, e.g. Newton-like Methods  Global Search, e.g. Genetic Algorithm

Forward Problem: Predicted Dispersion Curve

Initial Guess

31

Forward Modeling – State of the Art

 

2

1 1 1 1

T rr rz r rz zz z T T r z

r r r r r z z r z r r r z r

    

                                               u u u u D D D u D D D u u u D D D u I u

Discretize z direction Hankel Transform r direction

 

det  K

p

k c k    

Quadratic Eigenvalue Problem

Reducing the Problem Size: CFEM+PMDL

32

Complex-Length FEM (Finite Layers) Perfectly Matched discrete Layers (Halfspace)

Forward Modeling: CFEM vs. FEM

33

Error in dispersion curve

34

Inversion: Experimental Dispersion Curve

Experimental Dispersion Curve

Inverse Identification

240 Geophones 36 Geophones 12 Geophones

Experimental Dispersion Curve

1st (fundamental) Mode 2nd Mode 3rd Mode 4th Mode 5th Mode

35

Challenge

 No analytical derivative  Rough misfit function

 

2 experimantal predicted 1

Minimize

N i i i

E c c



 



 Existing approach: Finite Difference Method (FDM)  Expensive: Multiple computations of dispersion curve  Slow convergence: Oscillatory gradient

36

 Analytical Derivative

Proposed Derivative for Experimental Curve

L R L R

( ) ( )

j i j i

m k m k         K K    

† †

 Alternative Approach for Eigenvector

1 R 

 K P 

experimental

k k 

6

( 10 )

i

k k









 

37

 Analytical vs. FDM Gradient

Proposed Derivative

1 1 eff. 1 1

( )( / )( ) ( )( / )( )

j i j i

m k m k

   

        P K K K P P K K K P

† †

 FDM Gradient

Inversion Results: Synthetic Examples

38

 Analytical Gradient

Iterations (Existing) Iterations (Proposed) CPU Time (Proposed) CPU Time (Existing) 14 8 11.3 s 2884.6 s

Inversion Results: 14-Layer Soil Profile†

39 Analytical Gradient FDM Gradient

† Experimental data from: J Xia et al., J. Environ. Eng. Geophys., 5.3, 1-13 (2000)

Conclusions

 Discretization that perfectly preserves the impedance is possible

 Linear FEM with midpoint integration preserves impedance  Related to Crank-Nicolson discretization of the propagator  Preserves the evanescence in PML region

 Absorbing Boundary Conds.: Perfectly Matched Discrete Layers (PMDL)

 Exponential convergence  Link to other ABCs – we can get the best of both worlds  Facilitates stable ABCs for some backward propagating waves  Formally extensible to discrete periodic media  Open questions: Parameters of discretization for stability and accuracy

40

Conclusions (Contd.)

 Two-sided DtN Map

 Exponential convergence on the edges is possible with linear interpolation: Complex-length Finite Element Method (CFEM)  Impedance preserving discretization is the key!  Currently based on Padé approximant; could be further optimized  Open questions: further theoretical understanding; extensions to variable coefficients and higher dimensions?

 Guided Wave Inversion

 Forward modeling: a good application of CFEM  Approximate differentiation of the effective dispersion curve facilitates faster convergence and efficient gradient computation  Future work: Bayesian and hybrid inversion

41

This work is supported by National Science Foundation DMS-1016514, CMMI-1635291

Thank you!

Impedance Preserving Discretization

G (2006), Arbitrarily wide angle wave equations for complex media, CMAME

Perfectly Matched Discrete Layers

G, Tassoulas (2000), Continued fraction absorbing boundary conditions for the wave equation, J Comp. Acoustics. Asvadurov, Druskin, G, Knizhnerman (2003), On optimal finite-difference approximation of PML, SINUM. G, Lim (2006), Continued fraction absorbing boundary conditions for convex polygonal domains, IJNME. Thirunavukkarasu, G (2011), Absorbing boundary conditions for time-harmonic wave propagation in discretized domains, CMAME. Savadatti, G (2012), Accurate absorbing boundary conditions for anisotropic elastic media. parts 1/2, JCP.

Complex-length FEM

G, Druskin, Vaziri Astaneh (2016), Exponential Convergence through Linear Finite Element Discretization of Stratified Subdomains, JCP. Vaziri Astaneh, G (2016), Efficient Computation of Dispersion Curves for Multilayered Waveguides and Half-Spaces, CMAME.

Guided Wave Inversion

Vaziri Astaneh, G (2016), Improved Algorithms for Inversion of Surface Waves Using Multistation Analysis, Geoph. J. Intl.