[PPT] - Fast 3D Helmholtz Solvers for Seismic Inversion in the Frequency PowerPoint Presentation

SLIDE 1

Fast 3D Helmholtz Solvers for Seismic Inversion in the Frequency Domain

Russell J. Hewett

Mathematics & CMDA, Virginia Tech

Theory and Experience in Solving Inverse Problems in Geophysics Workshop Uppsala University

April 10, 2019

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 1 / 42

SLIDE 2

Collaborators

◮ Leonardo Zepeda-Nu˜

nez, Lawrence Berkeley National Lab

◮ Matthias Taus, TU Wien ◮ Laurent Demanet, MIT ◮ Adrien Scheuer, Universit`

e Catholique de Louvain

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 2 / 42

SLIDE 3

FWI in the Frequency Domain

PDE constrained optimization in frequency domain

◮ min J(m) = 1

2||d−F(m)||2 2 s.t. Lu = f

Advantages:

◮ No need to invert source time series

ˆ f(ω) = FFT(f(t))

◮ Only need specific frequency components RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 3 / 42

SLIDE 4

FWI in the Frequency Domain

PDE constrained optimization in frequency domain

◮ min J(m) = 1

2||d−F(m)||2 2 s.t. Lu = f

Advantages:

◮ Reduced memory and disk requirements in inverse problem

δm = − q, ∂ttu0T = − T q(x, t)∂ttu0(x, t)dt becomes δm = −

q, −ω2u0
Ω = −
ω

ˆ q(x, ω)−ω2ˆ u0(x, ω)

◮ Hybrid modeling: Use time-domain + DFT to achieve frequency

domain update

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 3 / 42

SLIDE 5

FWI in the Frequency Domain

PDE constrained optimization in frequency domain

◮ min J(m) = 1

2||d−F(m)||2 2 s.t. Lu = f

Advantages:

◮ Multiple simultaneous right-hand sides ◮ With a factorization based method, only need to Helmholtz operator

nce per domain

◮ Compare to explicit time-stepping: matvec required for each time step

for each source

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 3 / 42

SLIDE 6

FWI in the Frequency Domain

PDE constrained optimization in frequency domain

◮ min J(m) = 1

2||d−F(m)||2 2 s.t. Lu = f

Advantages:

◮ Heirarchichal frequency “sweeping” ⇒ Convergence guarantees (E. Beretta, M.V. de Hoop, F. Faucher, O. Scherzer (SIMA 2016)) RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 3 / 42

SLIDE 7

FWI in the Frequency Domain

100 200 300 400 500 20 40 60 80 100 120 140 100 200 300 400 500 20 40 60 80 100 120 140 100 200 300 400 500 20 40 60 80 100 120 140 Created with PySIT (www.pysit.org). RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 4 / 42

SLIDE 8

FWI in the Frequency Domain

PDE constrained optimization in frequency domain

◮ min J(m) = 1

2||d−F(m)||2 2 s.t. Lu = f

Challenges:

◮ Helmholtz in high frequency regime ◮ Helmholtz in 3D at high resolution ◮ Scalable Helmholtz in HPC environment RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 5 / 42

SLIDE 9

Motivation for Sweeping Solvers

Helmholtz at high frequency is hard Hu = (−ω2 − △)u = f + ABCs

◮ Frequency ω grows with n ◮ Computational load N scales with nd RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 6 / 42

SLIDE 10

Motivation for Sweeping Solvers

Helmholtz at high frequency is hard Hu = (−ω2 − △)u = f + ABCs

◮ Frequency ω grows with n ◮ Computational load N scales with nd

Classical dense direct methods in 3D

◮ memory-intensive ◮ hard to parallelize RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 6 / 42

SLIDE 11

Motivation for Sweeping Solvers

Helmholtz at high frequency is hard Hu = (−ω2 − △)u = f + ABCs

◮ Frequency ω grows with n ◮ Computational load N scales with nd

Classical dense direct methods in 3D

◮ memory-intensive ◮ hard to parallelize

Multigrid methods

◮ poor frequency scaling ◮ down-sampling oscillatory waves is hard RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 6 / 42

SLIDE 12

Motivation for Sweeping Solvers

Helmholtz at high frequency is hard Hu = (−ω2 − △)u = f + ABCs

◮ Frequency ω grows with n ◮ Computational load N scales with nd

Classical dense direct methods in 3D

◮ memory-intensive ◮ hard to parallelize

Multigrid methods

◮ poor frequency scaling ◮ down-sampling oscillatory waves is hard

Classical iterative schemes

◮ niter grows with ω RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 6 / 42

SLIDE 13

Sweeping Solvers and Domain Decomposition Methods

Sweeping Solvers/Preconditioners

◮ First O(N) claim (Engquist and Ying, 2010) ◮ First O(N) claim w/ domain decomposition (Stolk 2013) RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 7 / 42

SLIDE 14

Sweeping Solvers and Domain Decomposition Methods

Sweeping Solvers/Preconditioners

◮ First O(N) claim (Engquist and Ying, 2010) ◮ First O(N) claim w/ domain decomposition (Stolk 2013)

Other domain-decomposition methods (DDMs):

◮ Multifrontal w/ HSS compression (Xia, et al., 2013) ◮ Hierarchical Poincare-Steklov methods (Gillman, et al., 2014) ◮ Common challenges: ◮ Hazy scalability ◮ Issues with rough media RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 7 / 42

SLIDE 15

Sweeping Solvers and Domain Decomposition Methods

Sweeping Solvers/Preconditioners

◮ First O(N) claim (Engquist and Ying, 2010) ◮ First O(N) claim w/ domain decomposition (Stolk 2013)

Other domain-decomposition methods (DDMs):

◮ Multifrontal w/ HSS compression (Xia, et al., 2013) ◮ Hierarchical Poincare-Steklov methods (Gillman, et al., 2014) ◮ Common challenges: ◮ Hazy scalability ◮ Issues with rough media

Our approach: DDMs + sweeping w/ polarized traces

◮ Use direct methods distributed over tractable subproblems ◮ Glue with boundary integral formulations ◮ Embedded within iterative scheme RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 7 / 42

SLIDE 16

Take-home Messages

1. Polarized traces: domain decomposition done right

◮ Maximizes leveraging legacy direct solvers in the subdomains ◮ Maximal flexibility in parallel computation RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 8 / 42

SLIDE 17

Take-home Messages

1. Polarized traces: domain decomposition done right

◮ Maximizes leveraging legacy direct solvers in the subdomains ◮ Maximal flexibility in parallel computation

2. The number of iterations is

◮ Weakly dependent on the frequency ◮ Weakly dependent on the number L of subdomains ◮ Robust to roughness in background model RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 8 / 42

SLIDE 18

Take-home Messages

1. Polarized traces: domain decomposition done right

◮ Maximizes leveraging legacy direct solvers in the subdomains ◮ Maximal flexibility in parallel computation

2. The number of iterations is

◮ Weakly dependent on the frequency ◮ Weakly dependent on the number L of subdomains ◮ Robust to roughness in background model

3. Precondition interdomain communication with polarizing integral

conditions

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 8 / 42

SLIDE 19

Take-home Messages

1. Polarized traces: domain decomposition done right

◮ Maximizes leveraging legacy direct solvers in the subdomains ◮ Maximal flexibility in parallel computation

2. The number of iterations is

◮ Weakly dependent on the frequency ◮ Weakly dependent on the number L of subdomains ◮ Robust to roughness in background model

3. Precondition interdomain communication with polarizing integral

conditions

4. Can achieve sublinear complexity: better than O(RN)

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 8 / 42

SLIDE 20

Method Of Polarized Traces

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 9 / 42

SLIDE 21

Method Of Polarized Traces

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 9 / 42

SLIDE 22

Half-space Problem

f Γi,i+1 u

Polarization condition: 0 = −

Γ

G(x, y)∂nyu↑(y)dsy +

Γ

A System for Traces

◮ Assume local PDE is solved in the bulk ◮ Traces can be found by solving

M u = f =        v1

n

v2

1

v2

n

. . . vL

1

      

◮ M is constructed from dense Green’s function blocks. . . ◮ . . . non-trivial to invert

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 12 / 42

SLIDE 33

Polarization

◮ Annihilation relation:

◮ If u↑ is an up-going wavefield, then the annihilator relations are true on

the lower-half plane, i.e. G↓,ℓ

i

(u↑

−, u↑ 1) = 0,

for i ≥ 1.

◮ If u↓ is a down-going wavefield, then the annihilator relations are true

n the upper-half plane, i.e.

G↑,ℓ

i

(u↓

n, u↓ +) = 0,

for i ≤ nℓ

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 13 / 42

SLIDE 34

Polarization

1. Seek to solve M u = f

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 14 / 42

SLIDE 35

Polarization

1. Seek to solve M u = f
2. Transform to underdetermined system
M

M u↓ u↑

= −f

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 14 / 42

SLIDE 36

Polarization

1. Seek to solve M u = f
2. Transform to underdetermined system
M

M u↓ u↑

= −f
3. Constrain with annihilation relations

M M A↓ A↑ u↓ u↑

= −

f

RJH (Virginia Tech)

Polarized Traces Uppsala / April 10, 2019 14 / 42

SLIDE 37

Polarization

1. Seek to solve M u = f
2. Transform to underdetermined system
M

M u↓ u↑

= −f
3. Constrain with annihilation relations

M M A↓ A↑ u↓ u↑

= −

f

4. Additional minor transformations and permutations

D↓ U L D↑

u = f

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 14 / 42

SLIDE 38

Polarization

1. Seek to solve M u = f
2. Transform to underdetermined system
M

M u↓ u↑

= −f
3. Constrain with annihilation relations

M M A↓ A↑ u↓ u↑

= −

f

4. Additional minor transformations and permutations

D↓ U L D↑

u = f
5. Final system of equations

M u = f

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 14 / 42

SLIDE 39

Structure of Polarized SIE matrix M

M M

◮ . . . M is far easier to invert

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 15 / 42

SLIDE 40

Solving for Traces

◮ M =

D↓ D↑

+

U L

◮ Block diagonal, and block upper

and lower triangular

0.12
0.1
0.08
0.06
0.04
0.02

0.02 0.04 0.06 0.08

Iteration 2 (4 domain sweeps)

RJH (Virginia Tech) Polarized Traces Uppsala / April 10, 2019 17 / 42

SLIDE 44