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.
Introduction¶
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
where t is the time and \omega is the center frequency.
To see the connection with a Gaussian, it helps to define
which gives
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 = 1e3 center_frequency = 14.5 sigma_2 = 1 / (np.pi * center_frequency) ** 2 time = np.linspace(wavelet_width_in_seconds, 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.ylabel('amplitude') plt.title('Custom ricker wavelet') plt.show()
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 # Pwave velocity in m/s. vs = 1847.5 # Swave 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, max_frequency=max_frequency).create_mesh()
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 sideset receivers, see this tutorial.
# Receivers. nr = 5 # Number of receivers. depth = 200.0 # Depth of the receivers rx0 = 1010.0 # xvalue of first receiver. rx1 = 1410.0 # xvalue of last receiver. receivers = [ sc.receiver.cartesian.SideSetPoint2D( 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 2D, 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 dimensiondependent weights, where the weight in xdirection 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, source_time_function=stf)
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) sim.add_boundary_conditions(boundaries) sim.validate() sim
<salvus_flow.simple_config.simulation.Waveform object at 0x7f15e457f940>
Compare results¶
salvus_flow.api.run( site_name="local", input_file=sim, ranks=4, output_folder="outputcustom", 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 `outputcustom`. * 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 builtin Ricker wavelet.
source = sc.source.cartesian.VectorPoint2D( x=sx, y=sy, fx=fx, fy=fy, source_time_function=sc.source.stf.Ricker(center_frequency=center_frequency) ) sim.physics.wave_equation.point_source = [source] salvus_flow.api.run( site_name="local", input_file=sim, ranks=4, output_folder="outputbuiltin", 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 `outputbuiltin`. * 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 ['outputcustom', 'outputbuiltin']: 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.