# 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.

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)

``````

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)

``````
``````

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

Your browser does not support the video tag.

``````

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

``````
``````

Out[19]:

``````

It will contain the whole simulation history, just like it did before.

``````

In [20]:

``````
``````

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

``````
``````

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

``````
``````

In [22]:

``````
``````

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)

``````
``````

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

``````

Resuming the simulation should eventually trigger the second source.

``````

In [46]:

``````
``````

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

``````
``````

In [47]:

``````
``````

Out[47]:

``````

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

``````

In [48]:

``````
``````

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

``````
``````

In [49]:

``````
``````

Out[49]:

``````
``````

In [52]:

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)

``````
``````

``````