simpegSeis/notebooks/SeismicRefraction_model3.ipynb

In [1]:

# add simpegseisc directory to path
import os
if '..' not in os.sys.path:
    os.sys.path.append('..')

In [2]:

from simpegseis import *
%pylab inline
from SimPEG.Utils import mkvc

Populating the interactive namespace from numpy and matplotlib

WARNING: pylab import has clobbered these variables: ['linalg']
`%matplotlib` prevents importing * from pylab and numpy

In [11]:

# Step1: set time
dt = 0.8e-4
ntime = 1700
time = np.arange(ntime)*dt
# Step4: set mesh
cs = 0.5
cs = 0.5
hx = np.ones(500.)*cs
hy = np.ones(200.)*cs
mesh = Mesh.TensorMesh([hx, hy], 'CC')
# Step2: set Tx and Rx
tlag = 0.008
options={'tlag':tlag, 'fmain':100.} # You need to set waveform to set Tx
nrx = 99
xtemp_rx = np.arange(nrx)*2.5+2.5+mesh.vectorCCx.min()
xy_rx = np.c_[xtemp_rx, np.ones(nrx)*mesh.vectorCCy.max()]
rx = AcousticRx(xy_rx)

tx = AcousticSrc(np.r_[mesh.vectorCCx.min(), mesh.vectorCCy.max()], time, [rx], **options)
# Step3: set survey (pass txlist)
survey = SurveyAcoustic([tx])
wave = tx.RickerWavelet()
# Step5: set problem (pass mesh) and pair with survey
prob = AcousticProblemPML(mesh)
prob.pair(survey)
# Step6: set boundary
prob.setPMLBC(20, dt, bcflag='left', const=1.)
prob.storefield = True
# Step7: set velocity model and check stability
v = np.ones(mesh.nC)*1500.
v[mesh.gridCC[:,1]<30.] = 800.
v[mesh.gridCC[:,1]<0.] = 2300.
prob.stabilitycheck(v, time, 100.)

You are good to go:)
>> Stability information
   dt: 8.00e-05 s
   Optimal dt: 1.09e-04 s
   Cell per wavelength (G): 1.60e+01
   Optimal G: 1.60e+01

In [12]:

mesh.plotImage(v)

Out [12]:

(<matplotlib.collections.QuadMesh at 0x7fc0bee43090>,)

In [13]:

plt.plot(time, wave)

Out [13]:

[<matplotlib.lines.Line2D at 0x7fc0bedaf850>]

In [14]:

# Step8: run forward
U = prob.fields(v)

>> Start Computing Acoustic Wave
>> dt: 8.00e-05 s
>> Optimal dt: 1.09e-04 s
>> Main frequency, fmain: 1.00e+02 Hz
>> Cell per wavelength (G): 1.60e+01
  Src at (-124.75, -124.75):    1/   1
>>Elapsed time: 4.37e+01 s

In [15]:

# Step9: project data
data = survey.projectFields(U)

In [16]:

print 50./1500.

0.0333333333333

In [17]:

icount = 1699
fig, ax = plt.subplots(1,1, figsize = (8, 5))
extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min()-mesh.vectorCCy.max(), mesh.vectorCCy.max()-mesh.vectorCCy.max()]
ax.imshow(np.flipud(U[0][:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'binary', extent=extent)  
ax.imshow(np.flipud(v.reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'jet', extent=extent, alpha = 0.2)  

Out [17]:

<matplotlib.image.AxesImage at 0x7fc0bed48810>

In [18]:

print data[0].shape

(99, 1700)

In [28]:

from JSAnimation import IPython_display
from matplotlib import animation
fig, ax = plt.subplots(1,1, figsize = (8, 8))
ax.set_xlabel('Easting (m)', fontsize = 16)
ax.set_ylabel('Depth (m)', fontsize = 16)
extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min()-mesh.vectorCCy.max(), mesh.vectorCCy.max()-mesh.vectorCCy.max()]
ax.set_xlim(extent[:2])
ax.set_ylim(extent[2:])


nskip = 79
def animate(i_id):
    icount = i_id*nskip
    frame = ax.imshow(np.flipud(U[0][:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'binary', extent=extent)    
    ax.imshow(np.flipud(v.reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'jet', extent=extent, alpha = 0.2)  
    return frame

anim = animation.FuncAnimation(fig, animate, frames=20, interval=40, blit=True)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-28-894a9449f93f> in <module>()
     21     return frame
     22 
---> 23 anim = animation.FuncAnimation(fig, animate, init_func=init, frames=20, interval=40, blit=True)

/home/isisilon/.virtualenvs/main/local/lib/python2.7/site-packages/matplotlib/animation.pyc in __init__(self, fig, func, frames, init_func, fargs, save_count, **kwargs)
   1163         self._save_seq = []
   1164 
-> 1165         TimedAnimation.__init__(self, fig, **kwargs)
   1166 
   1167         # Need to reset the saved seq, since right now it will contain data

/home/isisilon/.virtualenvs/main/local/lib/python2.7/site-packages/matplotlib/animation.pyc in __init__(self, fig, interval, repeat_delay, repeat, event_source, *args, **kwargs)
   1007 
   1008         Animation.__init__(self, fig, event_source=event_source,
-> 1009                            *args, **kwargs)
   1010 
   1011     def _step(self, *args):

/home/isisilon/.virtualenvs/main/local/lib/python2.7/site-packages/matplotlib/animation.pyc in __init__(self, fig, event_source, blit)
    634                                                       self._stop)
    635         if self._blit:
--> 636             self._setup_blit()
    637 
    638     def _start(self, *args):

/home/isisilon/.virtualenvs/main/local/lib/python2.7/site-packages/matplotlib/animation.pyc in _setup_blit(self)
    905         self._resize_id = self._fig.canvas.mpl_connect('resize_event',
    906                                                        self._handle_resize)
--> 907         self._post_draw(None, self._blit)
    908 
    909     def _handle_resize(self, *args):

/home/isisilon/.virtualenvs/main/local/lib/python2.7/site-packages/matplotlib/animation.pyc in _post_draw(self, framedata, blit)
    870             self._blit_draw(self._drawn_artists, self._blit_cache)
    871         else:
--> 872             self._fig.canvas.draw_idle()
    873 
    874     # The rest of the code in this class is to facilitate easy blitting

/home/isisilon/.virtualenvs/main/local/lib/python2.7/site-packages/matplotlib/backend_bases.pyc in draw_idle(self, *args, **kwargs)
   2024         if not self._is_idle_drawing:
   2025             with self._idle_draw_cntx():
-> 2026                 self.draw(*args, **kwargs)
   2027 
   2028     def draw_cursor(self, event):

/home/isisilon/.virtualenvs/main/local/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in draw(self)
    472 
    473         try:
--> 474             self.figure.draw(self.renderer)
    475         finally:
    476             RendererAgg.lock.release()

/home/isisilon/.virtualenvs/main/local/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     59     def draw_wrapper(artist, renderer, *args, **kwargs):
     60         before(artist, renderer)
---> 61         draw(artist, renderer, *args, **kwargs)
     62         after(artist, renderer)
     63 

/home/isisilon/.virtualenvs/main/local/lib/python2.7/site-packages/matplotlib/figure.pyc in draw(self, renderer)
   1163 
   1164         self._cachedRenderer = renderer
-> 1165         self.canvas.draw_event(renderer)
   1166 
   1167     def draw_artist(self, a):

/home/isisilon/.virtualenvs/main/local/lib/python2.7/site-packages/matplotlib/backend_bases.pyc in draw_event(self, renderer)
   1807         s = 'draw_event'
   1808         event = DrawEvent(s, self, renderer)
-> 1809         self.callbacks.process(s, event)
   1810 
   1811     def resize_event(self):

/home/isisilon/.virtualenvs/main/local/lib/python2.7/site-packages/matplotlib/cbook.pyc in process(self, s, *args, **kwargs)
    561             for cid, proxy in list(six.iteritems(self.callbacks[s])):
    562                 try:
--> 563                     proxy(*args, **kwargs)
    564                 except ReferenceError:
    565                     self._remove_proxy(proxy)

/home/isisilon/.virtualenvs/main/local/lib/python2.7/site-packages/matplotlib/cbook.pyc in __call__(self, *args, **kwargs)
    428             mtd = self.func
    429         # invoke the callable and return the result
--> 430         return mtd(*args, **kwargs)
    431 
    432     def __eq__(self, other):

/home/isisilon/.virtualenvs/main/local/lib/python2.7/site-packages/matplotlib/animation.pyc in _start(self, *args)
    646 
    647         # Now do any initial draw
--> 648         self._init_draw()
    649 
    650         # Add our callback for stepping the animation and

/home/isisilon/.virtualenvs/main/local/lib/python2.7/site-packages/matplotlib/animation.pyc in _init_draw(self)
   1194 
   1195         else:
-> 1196             self._drawn_artists = self._init_func()
   1197             if self._blit:
   1198                 for a in self._drawn_artists:

<ipython-input-28-894a9449f93f> in init()
     11 # initialization function: plot the background of each frame
     12 def init():
---> 13     line.set_data([], [])
     14     return (line,)
     15 

NameError: global name 'line' is not defined

In [40]:

from JSAnimation import IPython_display
from matplotlib import animation
from IPython.display import HTML
fig, ax = plt.subplots(1,1, figsize = (8, 8))
ax.set_xlabel('Easting (m)', fontsize = 16)
ax.set_ylabel('Depth (m)', fontsize = 16)
extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min()-mesh.vectorCCy.max(), mesh.vectorCCy.max()-mesh.vectorCCy.max()]
ax.set_xlim(extent[:2])
ax.set_ylim(extent[2:])


nskip = 79
def animate(i_id):
    icount = i_id*nskip
    frame = ax.imshow(np.flipud(U[0][:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'binary', extent=extent)    
    ax.imshow(np.flipud(v.reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'jet', extent=extent, alpha = 0.2)  
    return frame

animation.FuncAnimation(fig, animate, frames=20, interval=40, blit=False)

Out [40]:

– + Once Loop Reflect

In [ ]:

In [20]:

def traveltime(x, h, v):
    d = np.sqrt((0.5*x)**2+h**2)
    t = 2*(d/v)
    return t

In [21]:

print traveltime(2.5, 20., 800.)

0.0500975610684

In [22]:

def intercTime(x, h, v1, v2):
    out = 2*h*np.sqrt(v2**2-v1**2)/(v1*v2)

In [23]:

# print intercTime(0., 20., 800.)

In [24]:

print (2*20)/800.

0.05

In [25]:

extent = [0,2, 0, 2]
fig, ax = plt.subplots(1,1, figsize = (8, 5))
ax.imshow(data[0].T, cmap = 'RdBu', extent=extent)  

Out [25]:

<matplotlib.image.AxesImage at 0x7fc0a1a46050>

In [26]:

fig, ax = plt.subplots(1,1, figsize=(12, 4))
ax.plot(time+tlag, data[0][0,:])
ax.grid(True)

In [ ]:

In [ ]: