Transcript weilib: wave-equation imaging

weilib: wave-equation imaging
Paul Sava
[email protected]
My goal
• Quick introduction to the library
• Not much about the imaging science
• Ask the right questions
[email protected]
Library goals
• Flexible research
– Code recycling
– Efficiency vs. flexibility
• Cluster access
– Parallelization with MPI, OMP
• Large problems
– 3-D data
[email protected]
Characteristics
• Programming
– Fortran 90
– OO
• Problem types
– Shot-receiver migration
– Shot-profile migration
• Versions
[email protected]
Linear operators
Lm  d
m
L
d
call operator_init(…)
stat=operator(adj,add,m,d)
[email protected]
WE imaging
• Large datasets (m,d)
– Data partition
Lm  d
• Complicated operators (L)
– Operator partition
[email protected]
The problem
IMAGE (m)
w
DATA (d)
z
mLd
*
[email protected]
Cluster distribution
w
MPI
z
mi  L*d i
[email protected]
Outline
• Overview
• Objects
– Partitioning
– I/O
• Operators
– WE migration
– WE MVA
• Examples
[email protected]
WEI initialization
…
call
call
call
call
call
call
…
sep_init()
weipar_init()
weipar_report()
weitag_init()
weipro_init()
weipro_bcast()
[email protected]
WEI objects & procedures
• Objects
–
–
–
–
Data
Image
Slowness
Wavefield
• Procedures
–
–
–
–
file tag
dimensions
I/O
MPI communication
[email protected]
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(:,:,:,:,:)
[email protected]
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)
…
[email protected]
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
[email protected]
WEI workflow example
…
call
call
call
…
call
…
call
…
weiS_mpicopy(SS)
weiR_mpizero(RR)
weiD_mpisprd(DD)
weimigof(…)
weiR_mpisumm(RR)
[email protected]
I/O blocking
w
z
[email protected]
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
…
[email protected]
Shared memory parallelization
w
z
OMP
[email protected]
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
[email protected]
Outline
• Overview
• Objects
– Partitioning
– I/O
• Operators
– WE migration
– WE MVA
• Examples
[email protected]
WE migration
imaging
frequency
extrapolation
IMAGE (m)
DATA (d)
depth
[email protected]
Extrapolation
frequency
Wz
Wz+z
Wz  z  Wz e
dk z
k z  k z0 
ds
ikz z
s  s0 
s0
depth
[email protected]
Imaging
frequency
Wz  z 
 Rz  z
Wz+z
Rz+z
depth
[email protected]
WE migration
dk z
i
ds
ik z0 z
Wz  z  Wz e
e
Wz  z 
 Rz  z
 s  s0 z
s0
[email protected]
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
[email protected]
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
[email protected]
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
[email protected]
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
[email protected]
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
[email protected]
Outline
• Overview
• Objects
– Partitioning
– I/O
• Operators
– WE migration
– WE MVA
• Examples
[email protected]
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
[email protected]
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
[email protected]
Outline
• Overview
• Objects
– Partitioning
– I/O
• Operators
– WE migration
– WE MVA
• Examples
[email protected]
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
[email protected]
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
[email protected]
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)
[email protected]
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
[email protected]
Summary
•
•
•
•
Flexible, reusable f90 code
Cluster ready (MPI,OMP)
Standard operator interface
Applications
– WE datuming
– WE migration
– WE MVA
[email protected]
Resources ([email protected])
• Developers
– Paul, Bob
• Manual
– Marie
• Users
– Biondo, Antoine, Morgan, Brad, Nick, Alejandro
[email protected]