Basic Helmholtz solver

In [1]:

    

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
%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': 12,
}

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


    

In [3]:

    

cellSize    = 1             # m
freq        = 1e2           # Hz        (should this be multiplied by 2*pi?)
velocity    = 2500          # m/s
density     = 2700          # units of density
Q           = 2e2           # can be inf
nx          = 512           # count
nz          = 256           # count
freeSurf    = [True, False, False, False] # t r b l
dims        = (nz,nx)       # tuple
sLoc1       = [128,255]     # x, z
sLoc2       = [384,255]     # x, z
nPML        = 64
c           = np.ones(dims) * velocity + np.linspace(0,4000,nz).reshape((nz,1))
rho         = np.ones(dims) * density
systemConfig = {
    'dx':   cellSize,       # m
    'dz':   cellSize,       # m
    'freq':     freq,       # Hz
    'c':    flipud(c),      # m/s
    'rho':  flipud(rho),    # density
    'Q':        Q,          # can be inf
    'nx':       nx,         # count
    'nz':       nz,         # count
    'freeSurf': freeSurf,   # t r b l
    'nPML':     nPML,
}


    

In [4]:

    

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


    

In [5]:

    

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 [6]:

    

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


    

    




CPU times: user 3.8 s, sys: 151 ms, total: 3.95 s
Wall time: 252 ms

In [7]:

    

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

fig = plt.figure()

ax1 = plt.subplot(1,2,1, 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(1,2,2, 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')


    

    
Out[7]:





<matplotlib.text.Text at 0x7fb9d51b5c10>






    

In [8]:

    

fig = plt.figure()

ax3 = plt.subplot(1,1,1, 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')


    

    
Out[8]:





<matplotlib.text.Text at 0x7fb988551050>






    

In [8]: