This is a prototype of classes to run FD simulations of the wave equation for Fatiando a Terra.
Below each class implementation there are a few example use cases of what it can do. This was also an oportunity for me to play around with h5py and IPython's rich display.
Skip below to the example usage for some figures and animations.
All content can be freely used and adapted under the terms of the Creative Commons Attribution 4.0 International License.
In [1]:
%matplotlib inline
from __future__ import division
import cPickle as pickle
import numpy as np
from fatiando import utils
from fatiando.seismic.wavefd import Ricker, ElasticPSV, ElasticSH
from fatiando.vis import mpl
The sources classes represent source function wavelets with the given amplitude, central frequency, and time delay. They have a rich display feature for IPython notebooks that plots the source function. So just putting as the last element of a cell will plot the source from [-2T, 2T], where T is the period.
In [2]:
w = Ricker(amp=10, cf=20)
w
Out[2]:
Applying a delay shifts the wavelet in time:
In [3]:
w = Ricker(amp=-1, cf=50, delay=1)
w
Out[3]:
We can simulate horizontal S waves using the ElasticSH class. This class solves the 2D elastic wave equation.
First, I'll need to make the S wave velocity and density grids for the simulation.
In [4]:
shape = (300, 800)
velocity = 3000*np.ones(shape)
density = 2200*np.ones(shape)
Then, I can create my simulation class and add a Ricker wavelet point source to the middle of the grid.
In [5]:
sim = ElasticSH(velocity, density, spacing=10)
sim.add_point_source((shape[0]//2, shape[1]//2), Ricker(5, 20, 1/20))
The simulation is stored in an HDF5 file. It is created automatically on the current working directory and is given a unique name. For this run, the cache file is:
In [6]:
sim.cachefile.split('/')[-1]
Out[6]:
Running the simulation is as simple as calling run with the desired number of iterations. It will even print an ASCII progress bar. You can turn this off by setting verbose=False when making the simulation instance.
In [7]:
sim.run(200)
The simulation classes define some rich display for the simulations. By default, the display for the simulation object is the last time frame computed.
In [8]:
sim
Out[8]:
You can browse each time frame of the simulation using the explore method. This creates an IPython widget that allows you to interactively explore each time frame of the simulation using a slider. This will only work if you're running the notebook localy. In nbviewer, only the last frame should appear.
In [9]:
sim.explore()
Out[9]:
You can inspect individual frames by calling snapshot. If you are not using %matplotlib inline, you can pass embed=True to insert the image in the notebook instead of poping out a plot window.
In [10]:
sim.snapshot(frame=100)
You can get the computed displacement data for each frame of the simulation by indexing the simulation object. Indexing the simulation will fetch the simulated displacements from the cache file. You can think of the sim object as 3D numpy array. The first dimension is time, the second is z (depth) and the last is x.
To get the displacements for the 100th time step:
In [11]:
sim[100]
Out[11]:
As with numpy arrays, you can slice the simulation object:
In [12]:
sim[100, 120:180, 370:430]
Out[12]:
And so we can plot a slice of the simulation at a given time:
In [13]:
mpl.imshow(sim[100, 120:180, 370:430], cmap=mpl.cm.gray)
mpl.colorbar(pad=0)
Out[13]:
And we can also index the simulation over a time period to get a seismogram:
In [14]:
seismo = sim[:, 155, 395]
mpl.plot(np.arange(sim.size)*sim.dt, seismo, '-k')
Out[14]:
To animate the simulation, use the animate method. The code below will animate using matplotlib's animation package. It will show every 10 frames of the simulation. A plot window should popup with the animation if you are not using %matplotlib inline. Problem with this is that it won't show on the notebook so you can't put it on nbviewer.
anim = sim.animate(every=10, cutoff=1, interval=100)
In [15]:
sim.animate(every=5, cutoff=1, blit=True, fps=20, dpi=50, embed=True)
Out[15]:
To resume the simulation, call run again.
In [16]:
sim.run(150)
In [17]:
sim
Out[17]:
You can reload the simulation from the cache file at a later time. So as long as you have the cache file, your simulation is safe.
Let's store the cache file name and delete the simulation object.
In [18]:
fname = sim.cachefile
del sim
Now we can use the from_cache method to restore the simulation object.
In [19]:
reloaded = ElasticSH.from_cache(fname)
reloaded
Out[19]:
It will contain the whole simulation history, just like it did before.
In [20]:
reloaded.snapshot(frame=100, embed=True)
Out[20]:
And the best thing is, you can even resume the simulation as if nothing had happened!
WARNING: this will write to the cache file so the old sim object might not work anymore.
In [21]:
reloaded.run(200)
In [22]:
reloaded
Out[22]:
The cache file stores the source functions as well. So even if you interupt the simulation before a source goes off, you can be sure it will be triggered once you reload the simulation from cache.
In [41]:
sim = ElasticSH(velocity, density, spacing=10)
sim.add_point_source((shape[0]//2, shape[1]//2), Ricker(5, 20, 1/20))
# Add a delayed source
sim.add_point_source((shape[0]//3, shape[1]//3), Ricker(5, 20, 0.3))
In [42]:
sim.run(200)
Check that the second source has not been fired:
In [43]:
sim
Out[43]:
Now we delete the simulation and reload it from the cache.
In [44]:
fname = sim.cachefile
del sim
In [45]:
reloaded = ElasticSH.from_cache(fname)
Resuming the simulation should eventually trigger the second source.
In [46]:
reloaded.run(200)
In [47]:
reloaded
Out[47]:
Lets have a look at a seismic section for this class:
In [48]:
reloaded.run(200)
In [49]:
reloaded
Out[49]:
In [52]:
section = reloaded[400:, 0, :]
mpl.figure(figsize=(10, 4))
mpl.imshow(section, origin='upper', cmap=mpl.cm.gray)
Out[52]:
In [57]:
shape = (300, 500)
pvel = 4000*np.ones(shape)
svel = 3000*np.ones(shape)
density = 2200*np.ones(shape)
The source function for this kind of wave needs to have an x and z component. When adding a point source, you only have to specify the dip of the source (with respect to the horizontal) and the source wavelet. The class will take care of projecting the source onto the x and z directions.
In [58]:
sim = ElasticPSV(pvel, svel, density, spacing=10)
sim.add_point_source((shape[0]//2, shape[1]//2), dip=45, wavelet=Ricker(5, 20, 1/20))
Running the simulation is the same as it was with ElasticSH.
In [59]:
sim.run(300)
The rich display features all work for PSV as well. However, PSV simulations would have 2 panels to show. It would be better to condense them into a single representation. A very convinient way to show P and S waves is to plot the divergence plus the curl of the displacement vector field. The curl will show the S waves and the divergence the P waves.
In [60]:
sim
Out[60]:
You can also use the interactive widget by calling explore and inspect each time step of the simulation with snapshot.
In [62]:
sim.snapshot(150)
Plotting methods in ElasticPSV can take some extra arguments. For example, you can use plottype=['vectors'] to plot the displacement vector field instead of the divergence + curl (i.e., plottype=['wavefield']). scale controls the exaggeration of the vector amplitudes.
In [63]:
sim.snapshot(250, plottype=['vectors'], scale=200)
The vector field is a good way of showing the difference in the vibration direction of P and S waves. Notice that the outer ring has displacement along the propagation direction (P waves), while in the inner ring it is perpendicular (S wave).
You can also combine different plot types.
In [64]:
sim.snapshot(250, plottype=['vectors', 'wavefield'], scale=200)
The same arguments will also work for explore and animate.
In [65]:
sim.animate(every=10, plottype=['vectors', 'wavefield'], scale=100, cutoff=0.1,
blit=True, embed=True, fps=10, dpi=50)
Out[65]:
Reloading from the cache also works and also preserves your sources.
In [75]:
sim = ElasticPSV(pvel, svel, density, spacing=10)
sim.add_point_source((shape[0]//2, shape[1]//2), dip=45, wavelet=Ricker(5, 20, 1/20))
# Add a delayed source deeper down
sim.add_point_source((2*shape[0]//3, shape[1]//3), dip=-45, wavelet=Ricker(5, 20, 0.25))
In [76]:
sim.run(200)
In [77]:
sim
Out[77]:
Reloading from cache and resuming will eventualy trigger the second source.
In [78]:
fname = sim.cachefile
del sim
In [79]:
reloaded = ElasticPSV.from_cache(fname)
In [80]:
reloaded.run(350)
In [81]:
reloaded.snapshot(-1, plottype=['wavefield', 'vectors'], scale=200)
Indexing an ElasticPSV object will get you two arrays, one for the x-displacement and one for the z-displacement.
In [107]:
ux, uz = reloaded[-1, 100:200, 100:400]
fig, axes = mpl.subplots(2, 1, figsize=(8, 6))
axes[0].set_title('ux')
axes[0].imshow(ux, cmap=mpl.cm.gray, origin='upper')
axes[1].set_title('uz')
axes[1].imshow(uz, cmap=mpl.cm.gray, origin='upper')
Out[107]:
ElasticPSV also has the a blast source. This source explodes in all directions, creating only P waves.
In [82]:
sim = ElasticPSV(pvel, svel, density, spacing=10)
sim.add_blast_source((shape[0]//2, shape[1]//2), wavelet=Ricker(3, 20, 1/20))
In [83]:
sim.run(300)
In [92]:
sim
Out[92]:
In [93]:
sim.animate(every=10, cutoff=0.05, plottype=['wavefield', 'vectors'], scale=100,
blit=True, embed=True, fps=10, dpi=50)
Out[93]: