Basic Helmholtz solver

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
%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]:

    

sf = SEGYFile(os.path.join(os.environ['HOME'],'Desktop/line10.vp'))
c = sf[:].T


    

In [4]:

    

cellSize    = 50            # m
freq        = 1e1           # Hz        (should this be multiplied by 2*pi?)
density     = 2700          # units of density
Q           = 2e2           # can be inf
nx          = c.shape[1]    # count
nz          = c.shape[0]    # count
freeSurf    = [False, False, False, False] # t r b l
dims        = (nz,nx)       # tuple
sLoc1       = [12000,3000]  # x, z
sLoc2       = [20000,3000]  # x, z
nPML        = 15
systemConfig = {
    'dx':   cellSize,       # m
    'dz':   cellSize,       # m
    'freq':     freq,       # Hz
    'c':    flipud(c),      # m/s
    'Q':        Q,          # can be inf
    'nx':       nx,         # count
    'nz':       nz,         # count
    'freeSurf': freeSurf,   # t r b l
    'nPML':     nPML,
}


    

In [5]:

    

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


    

In [6]:

    

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


    

In [7]:

    

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


    

    




CPU times: user 1.81 s, sys: 62 ms, total: 1.88 s
Wall time: 119 ms

In [8]:

    

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

fig = plt.figure()

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

ax1 = plt.subplot(4,1,2, aspect=1)
plt.set_cmap(cm.bwr)
im1 = mesh.plotImage(u1, vType='N', ax=ax1)
im1[0].set_clim(-amax, amax)
plt.title('Source 1')

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

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

fig.tight_layout()


    

In [8]: