next up previous [pdf]

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

\begin{displaymath}
{\frac{\partial^2 P}{\partial t^2}} =
{V^2(\mathbf{x}) {\...
...partial^2 P}{\partial x_2^2}\right) +
{F(\mathbf{x},t)}}\;,
\end{displaymath} (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:
\begin{displaymath}
\mathbf{P}_{t+\Delta t} = \left[ V^2(\mathbf{x}) \nabla^2 \...
...] \Delta t^2 + 2 \mathbf{P}_{t} - \mathbf{P}_{t-\Delta t} \;,
\end{displaymath} (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

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

\begin{displaymath}
{\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)
\end{displaymath} (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.
wave
[pdf] [png] [scons]

  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
vp,vx
Figure 3.
Vertical velocity (a) and horizontal velocity (b) in the Hess VTI model.
[pdf] [pdf] [png] [png] [scons]

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()


next up previous [pdf]

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

2013-09-17