Testing Environment for Eurus

Shaun Hadden & Brendan Smithyman | October, 2015

Imports


In [1]:
import numpy as np
from eurus import Eurus

In [2]:
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib

%matplotlib inline

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

Helper functions


In [3]:
class Source(object):
    
    def __init__(self, sc):
                
        self._x, self._z = np.mgrid[
            0:sc['dx']*sc['nx']:sc['dx'],
            0:sc['dz']*sc['nz']:sc['dz']
        ]
    
    def __call__(self, x, z):
        
        dist = np.sqrt((self._x - x)**2 + (self._z - z)**2)
        srcterm = 1.*(dist == dist.min())
        nz, nx = self._x.shape
        
        return np.hstack([srcterm.T.ravel() / srcterm.sum(), np.zeros(nx*nz, dtype=np.complex128)])

Modelling setup


In [4]:
# Geometry parameters
dx          = 1.
dz          = 1.
nx          = 100
nz          = 200

# Bulk parameters
velocity    = 2000.     * np.ones((nz,nx))
density     = 1.        * np.ones((nz,nx))

# Anisotropy parameters
theta       = 0.        * np.ones((nz,nx))
epsilon     = 0.2       * np.ones((nz,nx))
delta       = 0.2       * np.ones((nz,nx))

# Other parameters
freq        = 200.
freeSurf    = [False, False, False, False]
nPML        = 10


# Pack values into systemConfig dictionary
systemConfig = {
    'nx':       nx,
    'nz':       nz,
    'dx':       dx,
    'dz':       dz,

    'c':        velocity,
    'rho':      density,
    
    'theta':    theta,
    'eps':      epsilon,
    'delta':    delta,

    'freq':     freq,
    'freeSurf': freeSurf,
    'nPML':     nPML,
    'cPML':     1e3,
}

Testing


In [5]:
Ainv = Eurus(systemConfig)
src = Source(systemConfig)

In [6]:
plt.spy(Ainv.A, markersize=0.1, marker=',')


Out[6]:
<matplotlib.lines.Line2D at 0x10745f450>

In [7]:
q = src(50,100)
u = Ainv*q

In [8]:
clip = 1e-1

plotopts = {
    'extent':   [0, nx*dx, nz*dz, 0],
    'cmap':     cm.bwr,
    'vmin':     -clip,
    'vmax':     clip,
}

fig = plt.figure()

ax = fig.add_subplot(1,1,1, aspect=1)
plt.imshow(u[:nx*nz].reshape((nz,nx)).real, **plotopts)
plt.title('Real Wavefield')
plt.xlabel('X')
plt.ylabel('Z')


Out[8]:
<matplotlib.text.Text at 0x107440b90>

In [ ]:


In [ ]: