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           = inf           # can be inf
nx          = 512           # count
nz          = 512           # count
freeSurf    = [True, False, False, False] # t r b l
dims        = (nx,nz)       # tuple
sLoc        = [128,128]     # x, z
nPML        = 64
c           = fliplr(np.ones(dims) * velocity + np.linspace(0,4000,nz).reshape((1,nz)))
rho         = fliplr(np.ones(dims) * density)
systemConfig = {
    'dx':   cellSize,       # m
    'dz':   cellSize,       # m
    'freq':     freq,       # Hz
    'c':        c.T,          # m/s
    'rho':      rho.T,        # density
    'Q':        Q,          # can be inf
    'nx':       nx,         # count
    'nz':       nz,         # count
    'sLoc':     sLoc,       # x, z
    'freeSurf': freeSurf,   # t r b l
    'nPML':     nPML,
}


    

In [4]:

    

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


    

In [5]:

    

q       = np.zeros(mesh.nN)
qI      = Utils.closestPoints(mesh, sLoc, gridLoc='N')
q[qI]   = 1


    

In [6]:

    

%%time
u       = Ainv * q


    

    




CPU times: user 694 ms, sys: 26.7 ms, total: 721 ms
Wall time: 189 ms

In [7]:

    

amax    = abs(u).max()

ax1 = plt.subplot(1,2,1, aspect=1)
plt.set_cmap(cm.gray_r)
im1 = mesh.plotImage(q, vType='N', ax=ax1)
plt.title('Source location')

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


    

    
Out[7]:





<matplotlib.text.Text at 0x10c7355d0>






    

In [8]:

    

%%time
sLoc    = [128,256] # x, z
q       = np.zeros(mesh.nN)
qI      = Utils.closestPoints(mesh, sLoc, gridLoc='N')
q[qI]   = 1
u       = Ainv * q


    

    




CPU times: user 751 ms, sys: 18.3 ms, total: 770 ms
Wall time: 204 ms

In [9]:

    

amax    = abs(u).max()

ax1 = plt.subplot(1,2,1, aspect=1)
plt.set_cmap(cm.gray_r)
im1 = mesh.plotImage(q, vType='N', ax=ax1)
plt.title('Source location')

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


    

    
Out[9]:





<matplotlib.text.Text at 0x1ae9b0350>






    

In [10]:

    

%%time
sLoc    = [256,256] # x, z
q       = np.zeros(mesh.nN)
qI      = Utils.closestPoints(mesh, sLoc, gridLoc='N')
q[qI]   = 1
u       = Ainv * q


    

    




CPU times: user 750 ms, sys: 16.4 ms, total: 767 ms
Wall time: 203 ms

In [11]:

    

amax    = abs(u).max()

ax1 = plt.subplot(1,2,1, aspect=1)
plt.set_cmap(cm.gray_r)
im1 = mesh.plotImage(q, vType='N', ax=ax1)
plt.title('Source location')

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


    

    
Out[11]:





<matplotlib.text.Text at 0x1daf2c8d0>






    

In [11]: