Homework 2

Next: Programming part 1 Up: Homework 2 Previous: Theoretical part

# Computational part 1

In the first part of the computational assignment, you will investigate the numerical accuracy of a finite-difference eikonal solver.

 ```from rsf.proj import * from rsf.prog import RSFROOT import math # Program compilation ##################### proj = Project() # To do the coding assignment in Fortran, # comment the next line and uncomment the lines below prog = proj.Program('analytical.c') #prog = proj.Program('analytical.f90', # F90PATH=os.path.join(RSFROOT,'include'), # LIBS=['rsff90']+proj.get('LIBS')) exe = str(prog[0]) # Velocity model ######################################## v0 = 0.5 # velocity at the origin gz = 0.8 # vertical gradient gx = 0.6 # horizontal gradient nz = 601 # depth samples nx = 1201 # lateral samples dz = 0.5/(nz-1) # depth sampling dx = 1.0/(nx-1) # lateral sampling # v(z,x) = v0*/sqrt(1 - 2*v0*v0*(gz*z + gx*x)) ############################################## # !!! COMMENT BELOW !!! Flow('vel',None, ''' math n1=%d d1=%g n2=%d d2=%g output="%g/sqrt(1-%g*(%g*x1+%g*x2))" ''' % (nz,dz,nx,dx,v0,2*v0*v0,gz,gx)) # v(z,x) = v0 + gz*z + gx*x ########################### # !!! UNCOMMENT BELOW !!! #Flow('vel',None, # ''' # math n1=%d d1=%g n2=%d d2=%g # output="%g+%g*x1+%g*x2" # ''' % (nz,dz,nx,dx,v0,gz,gx)) Plot('vel', ''' grey color=j allpos=y screenht=6 screenratio=%g bias=%g wanttitle=n barreverse=y label1=Depth unit1=km label2=Lateral unit2=km scalebar=y barlabel=Velocity barunit="km/s" ''' % ((nz-1.)/(nx-1.),v0)) # Analytical traveltime ######################## # !!! CHANGE type=s to case=v BELOW !!! Flow('time',['vel',exe], ''' ./\${SOURCES[1]} v0=%g g1=%g g2=%g case=s ''' % (1./(v0*v0),-gz,-gx)) Plot('time', ''' contour scalebar=y plotcol=7 c0=0 dc=0.1 nc=30 screenht=6 screenratio=%g plotfat=5 wanttitle=n ''' % ((nz-1.)/(nx-1.))) # Plot wavefronts overlayed on the velocity model ################################################# Result('exact','vel time','Overlay') errs = [] spaces = [] for sample in range(1,nz//30): # Subsample velocity and analytical traveltime ############################################## vel = 'vel%d' % sample tim = 'tim%d' % sample window = 'window j1=%d j2=%d' % (sample,sample) Flow(vel,'vel',window) Flow(tim,'time',window) space = 'spc%d' % sample h = sample*math.hypot(dx,dz) Flow(space,None,'spike n1=1 mag=%g' % h) spaces.append(space) # Solve the eikonal equation ############################ eik = 'eik%d' % sample # !!!!!!!!! CHANGE BELOW !!!!!!!!!! Flow(eik,vel, ''' eikonal yshot=0 order=1 br1=%g br2=%g ''' % (20*dz,20*dx)) # Compute error ############### err = 'err%d' % sample Flow(err,[eik,tim], ''' math ref=\${SOURCES[1]} output="abs(input-ref)" | stack axis=2 norm=y | stack axis=1 norm=y ''') errs.append(err) # Concatenate results ##################### cat = 'cat axis=1 \${SOURCES[1:%d]}' % len(spaces) Flow('space',spaces,cat) Flow('logspace','space','math output="log(input)" ') Flow('err',errs,cat) Flow('logerr','err','math output="log(input)" ') graph = 'cmplx \${SOURCES[0:2]} | graph wanttitle=n ' for case in ('0','1'): Plot('err'+case,'space err',graph + ''' label1="Grid Spacing" unit1=km label2=Error unit2=s ''',stdin=0) Plot('logerr'+case,'logspace logerr',graph + ''' label1="Log(Grid Spacing)" unit1= label2="Log(Error)" unit2= grid=y ''',stdin=0) graph = graph + 'symbol=x plotcol=2 ' # Plot error versus spacing ########################### Plot('err','err0 err1','Overlay') # Plot error versus spacing on a logarithmic scale ################################################## Plot('logerr','logerr0 logerr1','Overlay') Result('err','err logerr','SideBySideAniso') End() ```

Figure 1 shows wavefronts in a medium with a constant gradient of slowness squared computed with an analytical two-point ray tracing formula from the class.

exact
Figure 1.
Wavefronts in a constant velocity gradient medium computed with an analytical formula.

By computing traveltimes numerically at different sampling intervals and comparing the numerical result with the analytical one, we can get an experimental estimate of the numerical error behavior. The error is shown in Figure 2. The right plot in Figure 2 displays the error in logarithmic coordinates. The slope of the line in these coordinates shows directly the rate of convergence of the numerical method. For example, a first-order accurate method should have a slope of one.

err
Figure 2.
Left: average error of the finite-difference eikonal solver as a function of grid spacing. Right: the same on a log-log plot. The slope of the curve on the log-log plot indicates the order of the numerical accuracy.

1. Change directory
```> cd hw2/eikonal
```
2. Run
```> scons view
```
to generate figures and display them on your screen.

3. In the SConstruct file, find the parameter that defines the order of accuracy for the eikonal solver. Change the order from to and recompute the results. Does the numerical accuracy change? What is the experimental order of accuracy?

 Homework 2

Next: Programming part 1 Up: Homework 2 Previous: Theoretical part

2013-09-19