Skip to content

Custom Source Time Functions

Jupyter Notebook

This page is automatically generated from a Jupyter Notebook.

Note that everything here is fully automatically tested against - and thus guaranteed to work - only for the latest versions of all the Salvus packages. So please update if anything works differently on your machine.

Click here for a version that is compatbile with the 0.9.X releases.


Salvus comes with a few wavelets that can readily be used as source time functions. However, sometimes it is necessary to define a custom source time function. This tutorial shows how to do that.

To this end, we are recreating Lamb's problem with our own implementation of the Ricker wavelet.

%matplotlib inline
%config Completer.use_jedi = False
%autosave 0
Autosave disabled
# Standard Python packages
import h5py
import numpy as np
import matplotlib.pyplot as plt

# Salvus packages
import salvus_flow.api
from salvus_mesh.simple_mesh import basic_mesh
import salvus_flow.simple_config as sc

Ricker wavelet as custom source time function

The Ricker wavelet is the second derivative of a Gaussian, and defined as

s(t) = \left(1 - 2\,t^2\,\pi^2\,\omega^2\right) \, \exp\left(-t^2\,\pi^2\,\omega^2\right),

where t is the time and \omega is the center frequency.

To see the connection with a Gaussian, it helps to define

\sigma = \left(\pi\,\omega\right)^{-1},

which gives

s(t) = \left(1 - \frac{2\,t^2}{\sigma^2}\right) \, \exp\left(-\frac{t^2}{\sigma^2}\right).

s(t) is centered around zero, so we either have to introduce a time shift, or start the simulation at a time t < 0.

There are two important things to notice:

  • To avoid artifacts in the wavefield the custom source function should always start smoothly from zero to be compatible with the homogeneous initial conditions.

  • You don't need to worry too much about the correct sampling rate, as the source time function will be resampled internally using the actual time step of the simulation. Just make sure that you have sufficiently many data points to avoid undersampling.

wavelet_width_in_seconds = 0.08
time_step_in_seconds = 1e-3
center_frequency = 14.5

sigma_2 = 1 / (np.pi * center_frequency) ** 2

time = np.linspace(-wavelet_width_in_seconds,
                   int((2 * wavelet_width_in_seconds / time_step_in_seconds)))

starttime_in_nanoseconds = 1e9 * time[0]
sampling_rate_in_hertz = 1.0 / time_step_in_seconds

wavelet = (1 - (2 * time ** 2) / sigma_2) * np.exp(-time ** 2 / sigma_2)

# plot the wavelet
plt.plot(time, wavelet)
plt.xlabel('time [s]')
plt.title('Custom ricker wavelet')


Simulation setup

To run a simulation with our custom source time function, we copy and paste all cells to setup the computational domain and mesh from here.

# Domain setup (m).
max_x = 2000.0       # Distance in meters.
max_y = 1000.0       # Distance in meters.
max_frequency = 25.0 # Frequency in Hz.

# Material properties.
vp  = 3200.0         # P-wave velocity in m/s.
vs  = 1847.5         # S-wave velocity in m/s.
rho = 2200.0         # Density in kg / m^3.
# Generate a mesh using the "simplified" meshing interface.
mesh = basic_mesh.CartesianHomogeneousIsotropicElastic2D(
    vp=vp, vs=vs, rho=rho, x_max=max_x, y_max=max_y,

We use the same receiver locations as in the previous tutorial, however, this time we define them relative to the surface using SideSetPoint2D receivers. For more information on side-set receivers, see this tutorial.

# Receivers.
nr = 5         # Number of receivers.
depth = -200.0 # Depth of the receivers
rx0 = 1010.0   # x-value of first receiver.
rx1 = 1410.0   # x-value of last receiver.
receivers = [
        point=np.array([x,max_y]), direction='y',
        side_set_name='y1', offset=depth,
        station_code=f"{_i:03d}", fields=["displacement"])
        for _i, x in enumerate(np.linspace(rx0, rx1, nr))]

Next, we create a point source at [1000.0, 500.0] that injects the Ricker wavelet we defined above. Because the simulation is in 2-D, the source time function needs to have 2 spatial components for x and y. Let's assume, we want to inject a force source in vertical direction, in this case we just multiply the wavelet with dimension-dependent weights, where the weight in x-direction is zero.

Next, we write the source time function to an hdf5 file, from which it will be read during the simulation. In addition to the array holding the source time function, we need to specify the start time in nano seconds(!) as well as the sampling rate in Hertz.

fx, fy = 0.0, -1e10 # Source components (Nm)
source_time_function = np.array([fx * wavelet, fy * wavelet]).T
with  h5py.File('stf.h5','w') as f:
    dset = f.create_dataset('stf', data=source_time_function)
    dset.attrs['sampling_rate'] = sampling_rate_in_hertz
    dset.attrs['starttime'] = starttime_in_nanoseconds

With the wavelet contained in the hdf5 file, we can create a Custom source time function object that points to this file. The directivity and amplitude is already included in the custom source time function, so the spatial weights of the VectorPoint2D source are set to 1.

# Sources.
sx, sy = 1000.0, 500.0 # Source position (m)

stf = sc.source.stf.Custom(dataset_name='/stf', filename='stf.h5')
custom_source = sc.source.cartesian.VectorPoint2D(
    x=sx, y=sy, fx=1.0, fy=1.0,

This is all we need to define the custom source time function. Of course, we can also inject multiple custom sources at arbitrary points. For this example, however, we are happy with a single point source, and assemble all information in the Waveform object to launch a simulation. To

sim = sc.simulation.Waveform(mesh=mesh, sources=custom_source, receivers=receivers)

sim.physics.wave_equation.end_time_in_seconds = 2.0

# write receivers to hdf5 format and only output every second time step
sim.output.point_data.format = "hdf5"
sim.output.point_data.sampling_interval_in_time_steps = 2

boundaries = sc.boundary.Absorbing(
    side_sets=["x0", "x1", "y0"],
    taper_amplitude=0.0, width_in_meters=0.0)



<salvus_flow.simple_config.simulation.Waveform object at 0x7f15e457f940>

Compare results
    site_name="local", input_file=sim, ranks=4,
    output_folder="output-custom", get_all=True, overwrite=True)
Job `job_1909101834390214_85eae723d2` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 0.10.3
  * Floating point size: 32

* Downloaded 51.4 KB of results to `output-custom`.
* Total run time: 1.70 seconds.
* Pure simulation time: 1.17 seconds.

<salvus_flow.sites.salvus_job.SalvusJob at 0x7f15e457d550>

Before we take a look at the results, let's change the source time function and rerun the same simulation with the built-in Ricker wavelet.

source = sc.source.cartesian.VectorPoint2D(
            x=sx, y=sy, fx=fx, fy=fy,

sim.physics.wave_equation.point_source = [source]
    site_name="local", input_file=sim, ranks=4,
    output_folder="output-built-in", get_all=True, overwrite=True)
Job `job_1909101834028226_00569cda78` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 0.10.3
  * Floating point size: 32

* Downloaded 51.8 KB of results to `output-built-in`.
* Total run time: 1.44 seconds.
* Pure simulation time: 1.19 seconds.

<salvus_flow.sites.salvus_job.SalvusJob at 0x7f15e4611978>

Now, we can compare the traces for both runs.

# Plot results.
fig, axes = plt.subplots(2, 5, sharex='col', sharey='row', figsize=(20,5))

for d in ['output-custom', 'output-built-in']:
    with h5py.File('{0}/receivers.h5'.format(d), 'r') as rec:

        coordinates = rec['coordinates_ELASTIC_point'][()]
        stations = rec['names_ELASTIC_point'][()]
        waveforms = rec['point']['displacement'][()]

        starttime = rec['point'].attrs['start_time_in_seconds']
        dt = 1/rec['point'].attrs['sampling_rate_in_hz']
        time = np.linspace(starttime, starttime+dt*waveforms.shape[2],waveforms.shape[2] )

        for r in range(waveforms.shape[0]):
            axes[0, r].plot(time, waveforms[r,0,:], label='cmp X, {0}'.format(stations[r]))
            axes[1, r].plot(time, waveforms[r,1,:], label='cmp Y, {0}'.format(stations[r]))
            axes[0, r].legend()
            axes[1, r].legend()

for a in axes[:, 0]:
    a.set_ylabel("Displacement (m)")
for a in axes[1, :]:
    a.set_xlabel("Time [s]")


As expected the traces look exactly the same.