Finite-difference wave simulation

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.

License

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 source functions

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

Elastic SH waves

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]:
'ElasticSH-8wqxVE.h5'

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)


|##################################################|100% Ran 200 iterations in 6.63163 seconds.

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]:
<function fatiando.seismic.wavefd.plot>

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]:
array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.]])

As with numpy arrays, you can slice the simulation object:


In [12]:
sim[100, 120:180, 370:430]


Out[12]:
array([[  2.60317298e-07,   8.74915863e-07,  -6.45496603e-07, ...,
         -8.18761962e-06,  -6.45496603e-07,   8.74915863e-07],
       [  8.74915863e-07,  -7.26305098e-07,  -8.73484007e-06, ...,
         -1.31911982e-05,  -8.73484007e-06,  -7.26305098e-07],
       [ -6.45496603e-07,  -8.73484007e-06,  -1.31072982e-05, ...,
          3.14237808e-05,  -1.31072982e-05,  -8.73484007e-06],
       ..., 
       [ -8.18761962e-06,  -1.31911982e-05,   3.14237808e-05, ...,
          1.59627536e-04,   3.14237808e-05,  -1.31911982e-05],
       [ -6.45496603e-07,  -8.73484007e-06,  -1.31072982e-05, ...,
          3.14237808e-05,  -1.31072982e-05,  -8.73484007e-06],
       [  8.74915863e-07,  -7.26305098e-07,  -8.73484007e-06, ...,
         -1.31911982e-05,  -8.73484007e-06,  -7.26305098e-07]])

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]:
<matplotlib.colorbar.Colorbar instance at 0x7f19262f6a70>

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]:
[<matplotlib.lines.Line2D at 0x7f1926776710>]

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)

For a more web-friendly version, use embed=True to insert a WebM video into the notebook. You'll need avconv or ffmpeg for this to work (the default is to use avconv; pass writer='ffmpeg' to use ffmpeg instead).


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)


|##################################################|100% Ran 150 iterations in 6.8539 seconds.

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)


|##################################################|100% Ran 200 iterations in 8.13035 seconds.

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)


|##################################################|100% Ran 200 iterations in 7.57739 seconds.

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)


|##################################################|100% Ran 200 iterations in 9.26461 seconds.

In [47]:
reloaded


Out[47]:

Lets have a look at a seismic section for this class:


In [48]:
reloaded.run(200)


|##################################################|100% Ran 200 iterations in 8.33862 seconds.

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]:
<matplotlib.image.AxesImage at 0x7f1925ea86d0>

Elastic P-SV waves

The class ElasticPSV simulates the coupled P and SV elastic waves and behaves pretty much the same as ElasticSH. The simulation is slower because it has to calculate the x and z displacements (ElasticSH only calculates y displacements).


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)


|##################################################|100% Ran 300 iterations in 12.5835 seconds.

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)