SEG Chevron 2014 Benchmark

Brendan Smithyman | September–October 2014

Imports and setup



In [1]:

    
import os
import numpy as np
import scipy.sparse as sp
from SimPEG import Utils
from zephyr import *
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from pygeo.segyread import SEGYFile
import scipy.interpolate as interpolate
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib



In [2]:

    
import matplotlib
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('png')
matplotlib.rcParams['savefig.dpi'] = 150 # Change this to adjust figure size

# Plotting options
font = {
    'family': 'Bitstream Vera Sans',
    'weight': 'normal',
    'size': 8,
}

matplotlib.rc('font', **font)



In [3]:

    
# Set the data directory here to point to your copy of the Chevron dataset
DATA_DIRECTORY = os.path.join(os.environ['HOME'],'Data/ChevronSEG2014')
SEGY_VERBOSE = True



In [4]:

    
# Functions to convert to local unit in km
lu = 'm' # label for figure axes

if lu == 'km':
    m2lu = lambda x: x * 1e-3
    km2lu = lambda x: m2lu(x) * 1e3
    xlabelspacing = 1
else:
    lu = 'm'
    m2lu = lambda x: x
    km2lu = lambda x: m2lu(x) * 1e3
    xlabelspacing = 2

# Function to make array
arr = lambda x: np.array(x)



In [5]:

    
# Load the velocity model
sfVelocity = SEGYFile(DATA_DIRECTORY + '/SEG14.Vpsmoothstarting.segy', endian='Big', verbose=SEGY_VERBOSE)

# Get an array containing the model
vel0 = sfVelocity[:].T

# Print out useful information
print('Model size: %d×%d'%(sfVelocity.ntr, sfVelocity.ns))
print('dz from headers: %4.3f %s'%(m2lu(sfVelocity.bhead['hdt'] * 1e-3),lu))









    



Detecting machine endianness...
Little.

Trying to create memory map...
Success. Using memory-mapped I/O.

Reading SEG-Y headers...
Read SEG-Y headers.
	3820 traces present.

Big endian specified... Not autodetecting.
Big endian != Little endian, therefore Foreign.
FORMAT == 1
             ...converting from IBM floating point.

Model size: 3820×1200
dz from headers: 5.000 m



In [6]:

    
# Grid size
DX, DZ = [12.5, 5.0] # m
DX_SRC = m2lu(25) # m
DX_REC = m2lu(25) # m

# Plotting extents
MODEL_EXTENT = km2lu(arr((0, 47.7375, 5.995, 0))) # x0, x1, z1, z0; all in km
WELL_EXTENT = (1500, 4500, m2lu(1000), m2lu(2450)) # v0, v1 in m/s; z0, z1 in km
WAVELET_EXTENT = (0, 0.250, -1.1, 1.1) # t0, t1 in ms; relative amps
WAVELET_SAMPR = 2. / 3e3 # sample rate in seconds
SPECTRUM_EXTENT = (0, 50, 0, 1.1) # f0, f1, p0, p1

# Array depths
SOURCE_DEPTH = m2lu(15) # m
REC_DEPTH = m2lu(15) # m

# Well position
WELL_X = m2lu(39375) # m
WELL_INDEX = 3150

# Velocity range
VMIN = 1500 # m/s
VMAX = 4500 # m/s

# Array parameters
SOURCES = 1600
CHANNELS = 321



In [7]:

    
# From: http://stackoverflow.com/a/312464
def chunks(l, n):
    for i in xrange(0, len(l), n):
        yield l[i:i+n]

Configuration



In [8]:

    
cellSize    = m2lu(20)      # m
freq        = 3e0           # Hz        (should this be multiplied by 2*pi?)
density     = 2700          # units of density
waterVel    = 1510          # m/s
waterRho    = 1000          # units of density
Q           = inf           # can be inf
freeSurf    = [True, False, False, False] # t r b l
sLoc1       = [24000,5950]  # x, z
sLoc2       = [32000,5950]  # x, z



In [9]:

    
Z, X        = np.mgrid[MODEL_EXTENT[3]:MODEL_EXTENT[2]+DZ:m2lu(DZ), MODEL_EXTENT[0]:MODEL_EXTENT[1]+DX:m2lu(DX)]
z, x        = np.mgrid[MODEL_EXTENT[3]:MODEL_EXTENT[2]:cellSize, MODEL_EXTENT[0]:MODEL_EXTENT[1]:cellSize]
nz, nx      = z.shape
c           = interpolate.griddata((Z.ravel(), X.ravel()), vel0.ravel(), (z.ravel(), x.ravel()), method='nearest').reshape((nz, nx))
waterLayer  = 1 * ((c-1510) < 1)
rho         = waterRho*waterLayer + density*(1-waterLayer)
Q           = waterLayer*1000 + (1 - waterLayer)*200
nPML        = np.int(c.max() / (cellSize * freq))



In [10]:

    
systemConfig = {
    'dx':   cellSize,       # m
    'dz':   cellSize,       # m
    'freq':     freq,       # Hz
    'c':    flipud(c),      # m/s
    'rho':      rho,        # kg/m^3
    'Q':        Q,          # can be inf
    'nx':       nx,         # count
    'nz':       nz,         # count
    'freeSurf': freeSurf,   # t r b l
    'nPML':     nPML,
}

Mesh + system initialization



In [11]:

    
#mesh, A = initHelmholtzSimPEG(systemConfig)
mesh, A = initHelmholtzNinePoint(systemConfig)
#Solver = SimPEG.SolverWrapI(sp.linalg.gmres, True, 1e0)
Solver = SimPEG.SolverWrapD(sp.linalg.splu)
Ainv = Solver(A)

Source initialization



In [12]:

    
q1      = np.zeros(mesh.nN)
qI1     = Utils.closestPoints(mesh, sLoc1, gridLoc='N')
q1[qI1]  = 1
q2      = np.zeros(mesh.nN)
qI2     = Utils.closestPoints(mesh, sLoc2, gridLoc='N')
q2[qI2]  = 1

Testing



In [13]:

    
%%time
u1 = Ainv * q1
u2 = Ainv * q2









    



CPU times: user 3.3 s, sys: 96.7 ms, total: 3.4 s
Wall time: 888 ms



In [14]:

    
amax = 5e-3*abs(u1).max()

fig = plt.figure()
ax0 = plt.subplot(1,1,1, aspect=1)
plt.set_cmap(cm.jet)
im0 = mesh.plotImage(fliplr(c.T), vType='CC', ax=ax0)
im0[0].set_clim(1500, 4200)
plt.title('Velocity')

fig = plt.figure()
ax1 = plt.subplot(1,1,1, aspect=1)
plt.set_cmap(cm.gray)
im1 = mesh.plotImage(fliplr(rho.T), vType='CC', ax=ax1)
im1[0].set_clim(waterRho, density)
plt.title('Density')

fig = plt.figure()
ax2 = plt.subplot(1,1,1, aspect=1)
plt.set_cmap(cm.bwr)
im2 = mesh.plotImage(u1, vType='N', ax=ax2)
im2[0].set_clim(-amax, amax)
plt.title('Source 1')

fig = plt.figure()
ax3 = plt.subplot(1,1,1, aspect=1)
plt.set_cmap(cm.bwr)
im3 = mesh.plotImage(u2, vType='N', ax=ax3)
im3[0].set_clim(-amax, amax)
plt.title('Source 2')

fig = plt.figure()
ax4 = plt.subplot(1,1,1, aspect=1)
plt.set_cmap(cm.bwr)
im4 = mesh.plotImage(u1*u2, vType='N', ax=ax4)
im4[0].set_clim(-1e-2 * amax**2, 1e-2 * amax**2)
plt.title('Sensitivity')









    Out[14]:





<matplotlib.text.Text at 0x104ddaa50>



In [15]: