Homework 1

Computational part

In this part, we will simulate and observe acoustic wave propagation in a simple velocity model shown in Figure 1a. A wave snapshot is shown in Figure 1b.

model,snap
Figure 1. (a) Velocity model for simple wave propagation experiments. (b) Wave snapshot with overlaid first-arrival wavefront.

1. Change directory to hw1/wave.
2. Run

   scons model.view
   

   to generate and view the velocity model from Figure 1a.
3. Run

   scons wave.vpl
   

   to generate and observe a propagating wave on your screen.
4. Run

   scons fronts.vpl
   

   to generate and observe a propagating first-arrival wavefront on your screen.
5. Run

   scons snap.view
   

   to generate and view a wave snapshot selected at 1 s as shown in Figure 1b.
6. Open the file SConstruct in your favorite editor. Your task is to make the following modifications in it:
   • Find a parameter responsible for selecting the time frame for the snapshot in Figure 1b. Modify it to increase the time from 1 s to your favorite point in the movie. Run scons snap.view again to verify your change.
   • Can you observe a geometrical part of the wave that is not captured by the first-arrival wavefront? What is its physical meaning?
   • Find a parameter in the SConstruct file responsible for the vertical smoothness of the model in Figure 1a. Modify it to increase the smoothness of the model in such a way that the first-arrival wavefront follows the wave geometry exactly. Run scons snap.view again to verify your change.
   • (EXTRA CREDIT) For extra credit, modify SConstruct to generate a movie of pictures like Figure 1b for a gradually increasing smoothness.

from rsf.proj import *

# Make a velocity model with a hyperbolic reflector
Flow('model',None,
     '''
     math n1=301 o1=-1 d1=0.01 output="sqrt(1+x1*x1)" |
     unif2 n1=201 d1=0.01 v00=1,2 |
     put label1=Depth unit1=km label2=Lateral unit2=km 
     label=Velocity unit=km/s | 
     smooth rect1=3
     ''')

# Plot model
Result('model',
       '''
       grey allpos=y title=Model bias=1
       scalebar=y barreverse=y 
       ''')
    
# Source wavelet
Flow('wavelet',None,
     '''
     spike nsp=1 n1=2000 d1=0.001 k1=201 |
     ricker1 frequency=10
     ''')

# Extended model (for absorbing boundaries)
Flow('left','model',
     '''
     window n2=1 | spray axis=2 n=50 o=-1.5 d=0.01 |
     math output="input*exp(-4*(-1-x2)^2)"
     ''')
Flow('right','model',
     '''
     window n2=1 f2=300 | spray axis=2 n=50 o=3.01 d=0.01 |
     math output="input*exp(-4*(3-x2)^2)"
     ''')
Flow('emodel2','left model right',
     'cat axis=2 ${SOURCES[1:3]}')

Flow('top','emodel2',
     '''
     window n1=1 | spray axis=1 n=50 o=-0.5 d=0.01 |
     math output="input*exp(-4*x1^2)"
     ''')
Flow('bottom','emodel2',
     '''
     window n1=1 f1=200 | spray axis=1 n=50 o=2.01 d=0.01 |
     math output="input*exp(-4*(2-x1)^2)"
     ''')
Flow('emodel','top emodel2 bottom',
     'cat axis=1 ${SOURCES[1:3]}')

# Source location
Flow('source',None,
     '''
     spike k1=101 k2=251
     n2=401 o2=-1.5 d2=0.01 label1=Depth   unit1=km
     n1=301 o1=-0.5 d1=0.01 label2=Lateral unit2=km
     ''')

# Modeling
exe = Program('wave.c')
Flow('wave','source %s wavelet emodel' % exe[0],
     '''
     ./${SOURCES[1]} wav=${SOURCES[2]} v=${SOURCES[3]} 
     jt=5 ft=200
     ''')

# Movie of wave snapshots
Plot('wave',
     '''
     window f1=50 f2=50 n1=201 n2=301 | 
     grey gainpanel=all title=Wave
     ''',view=1)

# Your favorite time moment
###########################
time = 1.0 # !!! MODIFY ME
###########################

# Wavefield snapshot
Plot('snap','wave',
     '''
     window f1=50 f2=50 n1=201 n2=301 n3=1 min3=%g |
     grey title="Wave Snapshot at %g s"
     label1=Depth unit1=km label2=Lateral unit2=km
     ''' % (time,time))

# First-arrival traveltime
Flow('first','model',
     'eikonal yshot=1 zshot=0.5 | add add=0.2')

# Movie of first-arrival wavefronts
fronts = []
for snap in range(180):
    front = 'front%d' % snap
    fronts.append(front)
    tsnap = 0.2+snap*0.01
    Plot(front,'first',
         'contour nc=1 c0=%g title="%g s" ' % (tsnap,tsnap))
Plot('fronts',fronts,'Movie',view=1)

# First-arrival wavefront
Plot('front','first',
     '''
     contour nc=1 c0=%g wanttitle=n wantaxis=n
     plotcol=3 plotfat=5 
     ''' % time)

# Overlay wavefront and traveltime
Result('snap','snap front','Overlay')

End()

Homework 1