Transcript weilib: wave-equation imaging

weilib: wave-equation imaging
Paul Sava
paul@sep.stanford.edu
My goal
• Quick introduction to the library
• Not much about the imaging science
• Ask the right questions
paul@sep.stanford.edu
Library goals
• Flexible research
– Code recycling
– Efficiency vs. flexibility
• Cluster access
– Parallelization with MPI, OMP
• Large problems
– 3-D data
paul@sep.stanford.edu
Characteristics
• Programming
– Fortran 90
– OO
• Problem types
– Shot-receiver migration
– Shot-profile migration
• Versions
paul@sep.stanford.edu
Linear operators
Lm  d
m
L
d
call operator_init(…)
stat=operator(adj,add,m,d)
paul@sep.stanford.edu
WE imaging
• Large datasets (m,d)
– Data partition
Lm  d
• Complicated operators (L)
– Operator partition
paul@sep.stanford.edu
The problem
IMAGE (m)
w
DATA (d)
z
mLd
*
paul@sep.stanford.edu
Cluster distribution
w
MPI
z
mi  L*d i
paul@sep.stanford.edu
Outline
• Overview
• Objects
– Partitioning
– I/O
• Operators
– WE migration
– WE MVA
• Examples
paul@sep.stanford.edu
WEI initialization
…
call
call
call
call
call
call
…
sep_init()
weipar_init()
weipar_report()
weitag_init()
weipro_init()
weipro_bcast()
paul@sep.stanford.edu
WEI objects & procedures
• Objects
–
–
–
–
Data
Image
Slowness
Wavefield
• Procedures
–
–
–
–
file tag
dimensions
I/O
MPI communication
paul@sep.stanford.edu
WEI data object
type wedata
character(len=128)
character(len=128)
integer
type(axis)
complex, pointer
end type wedata
::
::
::
::
::
tag
from
bites,esize
amx,amy,ahx,ahy,aw
bin(:,:,:,:,:)
paul@sep.stanford.edu
WEI data object
DD%tag=wtag%D
call weiD_init(DD)
call weiD_report(DD)
call weiD_allocate(DD)
call weiD_read(DD)
call weiD_write(DD)
call weiD_mpisprd(DD)
call weiD_mpigthr(DD)
…
paul@sep.stanford.edu
WEI migration workflow
data
image
slowness
data
Node 1
image
slowness
MPI
Node 2
MPI
data
data
slowness
image
slowness
data
slowness
image
Node 3
image
Node 4
paul@sep.stanford.edu
WEI workflow example
…
call
call
call
…
call
…
call
…
weiS_mpicopy(SS)
weiR_mpizero(RR)
weiD_mpisprd(DD)
weimigof(…)
weiR_mpisumm(RR)
paul@sep.stanford.edu
I/O blocking
w
z
paul@sep.stanford.edu
I/O blocking example
…
do iwb=1,nwb
call weiD_read(bW,iwb)
do izb=1,nzb
call weiR_read(RR,izb)
call weiS_read(SS,izb)
…
call weiR_write(RR,izb)
end do
end do
…
paul@sep.stanford.edu
Shared memory parallelization
w
z
OMP
paul@sep.stanford.edu
Shared memory parallelization
!$OMP PARALLEL
!$OMP DO &
!$OMP$ PRIVATE(iws,izs,ith)
do iws=1,nws
ith=omp_get_thread_num()+1
do izs=1,nzs
…
end do
end do
!$OMP END DO
!$OMP END PARALLEL
paul@sep.stanford.edu
Outline
• Overview
• Objects
– Partitioning
– I/O
• Operators
– WE migration
– WE MVA
• Examples
paul@sep.stanford.edu
WE migration
imaging
frequency
extrapolation
IMAGE (m)
DATA (d)
depth
paul@sep.stanford.edu
Extrapolation
frequency
Wz
Wz+z
Wz  z  Wz e
dk z
k z  k z0 
ds
ikz z
s  s0 
s0
depth
paul@sep.stanford.edu
Imaging
frequency
Wz  z 
 Rz  z
Wz+z
Rz+z
depth
paul@sep.stanford.edu
WE migration
dk z
i
ds
ik z0 z
Wz  z  Wz e
e
Wz  z 
 Rz  z
 s  s0 z
s0
paul@sep.stanford.edu
WCop: wavefield continuation operator
dk z
i
ds
ik z0 z
Wz  z  Wz e
e
Wz  z 
 Rz  z
 s  s0 z
s0
• Single reference slowness
• Multi reference slowness
– Wavefield interpolation
paul@sep.stanford.edu
SLop: slowness selection operator
dk z
i
ds
ik z0 z
Wz  z  Wz e
e
Wz  z 
 Rz  z
 s  s0 z
s0
• Min/Mean/Median
• Lloyd
paul@sep.stanford.edu
FKop: phase-shift operator
dk z
i
ds
ik z0 z
Wz  z  Wz e
e
Wz  z 
 Rz  z
 s  s0 z
s0
• Common-azimuth
• Narrow-azimuth
• Full prestack
paul@sep.stanford.edu
FXop: finite-difference operator
dk z
i
ds
ik z0 z
Wz  z  Wz e
e
Wz  z 
 Rz  z
 s  s0 z
s0
• Phase-shift
• Split-step Fourier
• Generalized Screen Propagator
paul@sep.stanford.edu
IGop: imaging operator
dk z
i
ds
ik z0 z
Wz  z  Wz e
e
Wz  z 
 Rz  z
 s  s0 z
s0
• Offset
• Ray parameter
paul@sep.stanford.edu
Outline
• Overview
• Objects
– Partitioning
– I/O
• Operators
– WE migration
– WE MVA
• Examples
paul@sep.stanford.edu
WE MVA
Wz  z  Wz e
dk z
i
ds
ik z0 z
e
s z
s0


dk
0 
z
Wz  z  Wz  z 1  i
sz 


ds
s0


dk z
W  W0 i
zs
ds s0
paul@sep.stanford.edu
SCop: scattering operator
Wz  z  Wz e
dk z
i
ds
ik z0 z
e
dk z
W  W0 i
ds
s z
s0
zs
s0
paul@sep.stanford.edu
Outline
• Overview
• Objects
– Partitioning
– I/O
• Operators
– WE migration
– WE MVA
• Examples
paul@sep.stanford.edu
WEI datuming (survey sinking)
call weidtm_init(
SLin=weiop_slo1_init,
WCin=weiop_mwc1_init,
FKin=weiop_wem3_init,
FXin=weiop_ssf3_init)
Slowness operator
stat = weidtm(adj,add, D, U,
SLop=weiop_slo1,
WCop=weiop_mwc1,
FKop=weiop_wem3,
FXop=weiop_ssf3)
Continuation operator
f-k operator
f-x operator
paul@sep.stanford.edu
WEI migration: example 1
call weimig_init(
SLin=weiop_slo1_init,
WCin=weiop_mwc1_init,
FKin=weiop_wem3_init,
FXin=weiop_ssf3_init,
IGin=weiop_pcig_init)
stat = weimig(adj,add, R, D,
SLop=weiop_slo1,
WCop=weiop_mwc1,
FKop=weiop_wem3,
FXop=weiop_ssf3,
IGop=weiop_pcig)
Imaging operator
paul@sep.stanford.edu
WEI migration: example 2
call weimig_init(
SLin=weiop_sllN_init,
WCin=weiop_mwcN_init,
FKin=weiop_cam3_init,
FXin=weiop_ssf3_init,
IGin=weiop_hcig_init)
stat = weimig(adj,add, R, D,
SLop=weiop_sllN,
WCop=weiop_mwcN,
FKop=weiop_cam3,
FXop=weiop_ssf3,
IGop=weiop_hcig)
paul@sep.stanford.edu
WEI MVA
call weimva_init(
SLin=weiop_sllN_init,
SCin=weiop_bor1_init,
WCin=weiop_mwcN_init,
FKin=weiop_wem3_init,
FXin=weiop_ssf3_init,
IGin=weiop_hcig_init)
stat = weimva(adj,add, dS,dR,
SLop=weiop_sllN,
SCop=weiop_bor1,
WCop=weiop_mwcN,
FKop=weiop_wem3,
FXop=weiop_ssf3,
IGop=weiop_hcig)
Scattering operator
paul@sep.stanford.edu
Summary
•
•
•
•
Flexible, reusable f90 code
Cluster ready (MPI,OMP)
Standard operator interface
Applications
– WE datuming
– WE migration
– WE MVA
paul@sep.stanford.edu
Resources (name@sep.stanford.edu)
• Developers
– Paul, Bob
• Manual
– Marie
• Users
– Biondo, Antoine, Morgan, Brad, Nick, Alejandro
paul@sep.stanford.edu