Homework 1

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

$${\frac{\partial^2 P}{\partial t^2}} = {V^2(\mathbf{x})\,{\... ...partial^2 P}{\partial x_2^2}\right) + {F(\mathbf{x},t)}}\;,$$ (15)

where $F(\mathbf{x},t)$ 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:

$$\mathbf{P}_{t+\Delta t} = \left[ V^2(\mathbf{x})\,\nabla^2 \... ...]\,\Delta t^2 + 2 \mathbf{P}_{t} - \mathbf{P}_{t-\Delta t} \;,$$ (16)

where $\mathbf{P}$ represents the propagating wavefield discretized at different time steps.

/* 2-D finite-difference acoustic wave propagation */
#include <rsf.h>

#ifdef _OPENMP
#include <omp.h>
#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<n2; i2++) {
	for (i1=0; i1<n1; i1++) {
	    u0[i2][i1]=0.0;
	    u1[i2][i1]=0.0;
	    u2[i2][i1]=0.0;
	    ud[i2][i1]=0.0;
	    vv[i2][i1] *= vv[i2][i1]*dt2;
	}
    }

    /* Time loop */
    for (it=0; it<nt; it++) {
	laplacian(u1,ud);

#ifdef _OPENMP
#pragma omp parallel for	         private(i2,i1)		         shared(ud,vv,ww,it,rr,u2,u1,u0)
#endif	
	for (i2=0; i2<n2; i2++) {
	    for (i1=0; i1<n1; i1++) {
		/* scale by velocity */
		ud[i2][i1] *= vv[i2][i1];
		/* inject wavelet */
		ud[i2][i1] += ww[it] * rr[i2][i1];
		/* time step */
		u2[i2][i1] = 
		    2*u1[i2][i1] 
		    - u0[i2][i1] 
		    + ud[i2][i1]; 
		u0[i2][i1] = u1[i2][i1];
		u1[i2][i1] = u2[i2][i1];
	    }
	}		
	
	/* write wavefield to output */
	if (it >= 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

#!/usr/bin/env python

import sys
import numpy
import m8r

class Laplacian:
    'Laplacian operator, 4th-order finite-difference'
    
    def __init__(self,d1,d2):
        self.c11 = 4.0*d1/3.0
        self.c12=  -d1/12.0
        self.c21 = 4.0*d2/3.0
        self.c22=  -d2/12.0
        self.c0  = -2.0*(self.c11+self.c12+self.c21+self.c22)
    
    def apply(self,uin,uout):
        n1,n2 = uin.shape
    
        uout[2:n1-2,2:n2-2] =           self.c11*(uin[1:n1-3,2:n2-2]+uin[3:n1-1,2:n2-2]) +           self.c12*(uin[0:n1-4,2:n2-2]+uin[4:n1  ,2:n2-2]) +           self.c21*(uin[2:n1-2,1:n2-3]+uin[2:n1-2,3:n2-1]) +           self.c22*(uin[2:n1-2,0:n2-4]+uin[2:n1-2,4:n2  ]) +           self.c0*uin[2:n1-2,2:n2-2]

par = m8r.Par()

# setup I/O files
Fr=m8r.Input()       # source position
Fo=m8r.Output()      # output wavefield

Fv=m8r.Input ("v")   # velocity
Fw=m8r.Input ("wav") # source wavefield


# Read/Write axes
a1 = Fr.axis(1); n1 = a1['n']; d1 = a1['d']
a2 = Fr.axis(2); n2 = a2['n']; d2 = a2['d']
at = Fw.axis(1); nt = at['n']; dt = at['d']

ft = par.int('ft',0)
jt = par.int('jt',0)

Fo.put('n3',(nt-ft)/jt)
Fo.put('d3',jt*dt)
Fo.put('o3',ft*dt)

dt2 =    dt*dt

# set Laplacian coefficients 
laplace = Laplacian(1.0/(d1*d1),1.0/(d2*d2))

# read wavelet, velocity & source position
ww = numpy.zeros(nt,'f');      Fw.read(ww)
vv = numpy.zeros([n2,n1],'f'); Fv.read(vv)
rr = numpy.zeros([n2,n1],'f'); Fr.read(rr)

# allocate temporary arrays
u0 = numpy.zeros([n2,n1],'f')
u1 = numpy.zeros([n2,n1],'f')
u2 = numpy.zeros([n2,n1],'f')
ud = numpy.zeros([n2,n1],'f')

vv = vv*vv*dt2

# Time loop
for it in range(nt):
    laplace.apply(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 == (it-ft)%jt:
        sys.stderr.write("%d" % it)
        Fo.write(u1)

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

$${\frac{\partial^2 P}{\partial t^2}} = V_1^2(\mathbf{x})\,... ...f{x})\,\frac{\partial^2 P}{\partial x_2^2} + F(\mathbf{x},t)$$ (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 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')

# UNCOMMENT BELOW IF YOU WANT TO USE FORTRAN
#exe = proj.Program('wave.f90',
#                   F90PATH=os.path.join(RSFROOT,'include'),
#                   LIBS=['rsff90']+proj.get('LIBS'))

# UNCOMMENT BELOW IF YOU WANT TO USE PYTHON
#exe = proj.Command('wave.exe','wave.py','cp $SOURCE $TARGET')
#AddPostAction(exe,Chmod(exe,0o755))


# 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='hessvti',
          server='https://s3.amazonaws.com',
          top='open.source.geoscience/open_data')
    
    # 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