Homework 1

Next: Completing the assignment Up: Homework 1 Previous: Computational part

# Programming part (extra credit)

For extra credit, you can modify the wave modeling program to include anisotropic wave propagation effects. The program below (slightly modified from the original version by Paul Sava) implements wave modeling with equation

 (15)

where is the source term. The implementation uses finite-difference discretization (second-order in time and fourth-order in space). Stepping in time involves the following computations:
 (16)

where represents the propagating wavefield discretized at different time steps.

 /* 2-D finite-difference acoustic wave propagation */ #include #ifdef _OPENMP #include #endif static int n1, n2; static float c0, c11, c21, c12, c22; static void laplacian(float **uin /* [n2][n1] */, float **uout /* [n2][n1] */) /* Laplacian operator, 4th-order finite-difference */ { int i1, i2; #ifdef _OPENMP #pragma omp parallel for private(i2,i1) shared(n2,n1,uout,uin,c11,c12,c21,c22,c0) #endif for (i2=2; i2 < n2-2; i2++) { for (i1=2; i1 < n1-2; i1++) { uout[i2][i1] = c11*(uin[i2][i1-1]+uin[i2][i1+1]) + c12*(uin[i2][i1-2]+uin[i2][i1+2]) + c21*(uin[i2-1][i1]+uin[i2+1][i1]) + c22*(uin[i2-2][i1]+uin[i2+2][i1]) + c0*uin[i2][i1]; } } } int main(int argc, char* argv[]) { int it,i1,i2; /* index variables */ int nt,n12,ft,jt; float dt,d1,d2,dt2; float *ww,**vv,**rr; float **u0,**u1,**u2,**ud; sf_file Fw,Fv,Fr,Fo; /* I/O files */ /* initialize Madagascar */ sf_init(argc,argv); /* initialize OpenMP support */ omp_init(); /* setup I/O files */ Fr = sf_input ("in"); /* source position */ Fo = sf_output("out"); /* output wavefield */ Fw = sf_input ("wav"); /* source wavelet */ Fv = sf_input ("v"); /* velocity */ /* Read/Write axes */ if (!sf_histint(Fr,"n1",&n1)) sf_error("No n1= in inp"); if (!sf_histint(Fr,"n2",&n2)) sf_error("No n2= in inp"); if (!sf_histfloat(Fr,"d1",&d1)) sf_error("No d1= in inp"); if (!sf_histfloat(Fr,"d2",&d2)) sf_error("No d2= in inp"); if (!sf_histint(Fw,"n1",&nt)) sf_error("No n1= in wav"); if (!sf_histfloat(Fw,"d1",&dt)) sf_error("No d1= in wav"); n12 = n1*n2; if (!sf_getint("ft",&ft)) ft=0; /* first recorded time */ if (!sf_getint("jt",&jt)) jt=1; /* time interval */ sf_putint(Fo,"n3",(nt-ft)/jt); sf_putfloat(Fo,"d3",jt*dt); sf_putfloat(Fo,"o3",ft*dt); dt2 = dt*dt; /* set Laplacian coefficients */ d1 = 1.0/(d1*d1); d2 = 1.0/(d2*d2); c11 = 4.0*d1/3.0; c12= -d1/12.0; c21 = 4.0*d2/3.0; c22= -d2/12.0; c0 = -2.0 * (c11+c12+c21+c22); /* read wavelet, velocity & source position */ ww=sf_floatalloc(nt); sf_floatread(ww ,nt ,Fw); vv=sf_floatalloc2(n1,n2); sf_floatread(vv[0],n12,Fv); rr=sf_floatalloc2(n1,n2); sf_floatread(rr[0],n12,Fr); /* allocate temporary arrays */ u0=sf_floatalloc2(n1,n2); u1=sf_floatalloc2(n1,n2); u2=sf_floatalloc2(n1,n2); ud=sf_floatalloc2(n1,n2); for (i2=0; i2= ft && 0 == (it-ft)%jt) { sf_warning("%d;",it+1); sf_floatwrite(u1[0],n12,Fo); } } sf_warning("."); exit (0); } 

 ! 2-D finite-difference acoustic wave propagation module laplace ! Laplacian operator, 4th-order finite-difference implicit none real :: c0, c11, c21, c12, c22 contains subroutine laplacian_set(d1,d2) real, intent (in) :: d1,d2 c11 = 4.0*d1/3.0 c12= -d1/12.0 c21 = 4.0*d2/3.0 c22= -d2/12.0 c0 = -2.0 * (c11+c12+c21+c22) end subroutine laplacian_set subroutine laplacian(uin, uout) real, dimension (:,:), intent (in) :: uin real, dimension (:,:), intent (out) :: uout integer n1, n2 n1 = size(uin,1) n2 = size(uin,2) uout(3:n1-2,3:n2-2) = & c11*(uin(2:n1-3,3:n2-2) + uin(4:n1-1,3:n2-2)) + & c12*(uin(1:n1-4,3:n2-2) + uin(5:n1, 3:n2-2)) + & c21*(uin(3:n1-2,2:n2-3) + uin(3:n1-2,4:n2-1)) + & c22*(uin(3:n1-2,1:n2-4) + uin(3:n1-2,5:n2 )) + & c0*uin(3:n1-2,3:n2-2) end subroutine laplacian end module laplace program Wave use rsf use laplace implicit none integer :: it, nt, ft, jt, n1, n2 real :: dt, d1, d2, dt2 real, dimension (:), allocatable :: ww real, dimension (:,:), allocatable :: vv, rr real, dimension (:,:), allocatable :: u0,u1,u2,ud type(file) :: Fw,Fv,Fr,Fo ! I/O files call sf_init() ! initialize Madagascar ! setup I/O files Fr = rsf_input ("in") ! source position Fo = rsf_output("out") ! output wavefield Fw = rsf_input ("wav") ! source wavelet Fv = rsf_input ("v") ! velocity ! Read/Write axes call from_par(Fr,"n1",n1) call from_par(Fr,"n2",n2) call from_par(Fr,"d1",d1) call from_par(Fr,"d2",d2) call from_par(Fw,"n1",nt) call from_par(Fw,"d1",dt) call from_par("ft",ft,0) ! first recorded time call from_par("jt",jt,0) ! time interval call to_par(Fo,"n3",(nt-ft)/jt) call to_par(Fo,"d3",jt*dt) call to_par(Fo,"o3",ft*dt) dt2 = dt*dt ft = ft+1 ! set Laplacian coefficients call laplacian_set(1.0/(d1*d1),1.0/(d2*d2)) ! read wavelet, velocity & source position allocate (ww(nt), vv(n1,n2), rr(n1,n2)) call rsf_read(Fw,ww) call rsf_read(Fv,vv) call rsf_read(Fr,rr) ! allocate temporary arrays */ allocate (u0(n1,n2), u1(n1,n2), u2(n1,n2), ud(n1,n2)) u0=0.0 u1=0.0 u2=0.0 ud=0.0 vv = vv*vv*dt2 ! Time loop do it=1,nt call laplacian(u1,ud) ! scale by velocity ud = ud*vv ! inject wavelet ud = ud + ww(it) * rr ! time step u2 = 2*u1 - u0 + ud u0 = u1 u1 = u2 ! write wavefield to output if (it >= ft .and. 0 == mod(it-ft,jt)) then write(0,'(a,i4)',advance='no') achar(13), it call rsf_write(Fo,u1) end if end do write (0,*) call exit(0) end program Wave 

Your task is to modify the code to implement elliptically-anisotropic wave propagation according to equation

 (17)

You can test your implementation using a constant velocity example shown in Figure 2.

wave
Figure 2.
Wavefield snapshot for propagation from a point-source in a homogeneous medium. Modify the code to make wave propagation anisotropic.

1. Change directory to geo384w/hw1/code
2. Run
scons wave.vpl

to compile and run the program and to observe a propagating wave on your screen.
3. Open the file wave.c in your favorite editor and modify it to implement the wave operator from equation (17).
4. Run
scons wave.vpl

again to compile and test your program.
5. (EXTRA EXTRA CREDIT) For an additional test, modify the file SConstruct and run appropriate commands to output a snapshot of wave propagation in the Hess VTI model, shown in Figure 3. Modify the file paper.tex to include your figure.

vp,vx
Figure 3.
Vertical velocity (a) and horizontal velocity (b) in the Hess VTI model.

 from rsf.proj import * from rsf.prog import RSFROOT # Program compilation ##################### proj = Project() # To do the coding assignment in Fortran, # comment the next line and uncomment the lines below #exe = proj.Program('wave.c') exe = proj.Program('wave.f90', F90PATH=os.path.join(RSFROOT,'include'), LIBS=['rsff90']+proj.get('LIBS')) # Constant velocity test ######################## # Source wavelet Flow('wavelet',None, ''' spike n1=1000 d1=0.001 k1=201 | ricker1 frequency=10 ''') # Source location Flow('source',None, ''' spike n1=201 n2=301 d1=0.01 d2=0.01 label1=x1 unit1=km label2=x2 unit2=km k1=101 k2=151 ''') # Velocity model Flow('v1','source','math output=1') Flow('v2','source','math output=1.5') # Modeling Flow('wave','source %s wavelet v1 v2' % exe[0], ''' ./${SOURCES[1]} wav=${SOURCES[2]} v=${SOURCES[3]} vx=${SOURCES[4]} ft=200 jt=5 ''') Plot('wave','grey gainpanel=all title=Wave',view=1) Result('wave', ''' window n3=1 min3=0.9 | grey title=Wave screenht=8 screenwd=12 ''') # Download Hess VTI model ######################### zcat = WhereIs('gzcat') or WhereIs('zcat') for case in ('vp','epsilon'): sgy = 'timodel_%s.segy' % case sgyz = sgy + '.gz' Fetch(sgyz,dir='Hess_VTI', server='ftp://software.seg.org', top='pub/datasets/2D') # Uncompress Flow(sgy,sgyz,zcat + ' $SOURCE',stdin=0) # Convert to RSF format Flow(case,sgy, ''' segyread read=data | window j1=2 j2=2 | put d1=40 d2=40 unit1=ft label1=Depth unit2=ft label2=Distance ''') # Horizontal velocity Flow('vx','vp epsilon', 'math e=${SOURCES[1]} output="input*sqrt(1+2*e)"') for case in ('vp','vx'): Result(case, ''' grey color=j pclip=100 allpos=y bias=5000 scalebar=y barreverse=y wanttitle=n barlabel=Velocity barunit=ft/s screenht=5 screenwd=12 labelsz=6 ''') Flow('hsource','vp','spike k1=300 k2=900') Flow('hess','hsource %s wavelet vp vx' % exe[0], ''' ./${SOURCES[1]} wav=${SOURCES[2]} v=${SOURCES[3]} vx=${SOURCES[4]} ft=200 jt=5 ''') End() 

 Homework 1

Next: Completing the assignment Up: Homework 1 Previous: Computational part

2013-09-17