Simple example
SeisCL is of program that performs forward and adjoint seismic modeling for imaging and inversion. In this first example, we show how to use SeisCL to perform basic seismic modeling with the python wrapper. This simple example will solve the 2D isotropic elastic wave equation.
The first thing to do is to load the SeisCL class and relevant modules:
[1]:
from SeisCL import SeisCL
import numpy as np
import matplotlib.pyplot as plt
import os
SeisCL is controlled by a single class in python that will group all relevant information to perform forward and adjoint modeling. The class creates all parameters with default values: they can be set at inialization or can be changed later. Note that all the available parameters are defined in the docstring of the SeisCL class.
[2]:
seis = SeisCL()
The class creates all parameters with default values: they can be set at inialization or can be changed later. Note that all the available parameters are defined in the docstring of the SeisCL class, which can be accessed by using the help function:
[3]:
help(seis.__init__)
Help on method __init__ in module SeisCL.SeisCL:
__init__(N: numpy.ndarray = None, ND: int = 2, dh: float = 10, dt: float = 0.0008, NT: int = 875, L: int = 0, f0: float = 15, FL: numpy.ndarray = array(15), FDORDER: int = 8, MAXRELERROR: int = 1, FP16: int = 0, src_pos_all: numpy.ndarray = array([], shape=(5, 0), dtype=float64), rec_pos_all: numpy.ndarray = array([], shape=(8, 0), dtype=float64), src_all: numpy.ndarray = None, freesurf: int = 0, abs_type: int = 1, VPPML: float = 3500, NPOWER: float = 2, FPML: float = 15, K_MAX_CPML: float = 2, nab: int = 16, abpc: float = 6, with_docker: bool = False, with_mpi: bool = False, NP: int = 1, pref_device_type: int = 4, MPI_NPROC_SHOT: int = 1, nmax_dev: int = 1, no_use_GPUs: numpy.ndarray = array([-1]), gradout: int = 0, Hout: int = 0, gradsrcout: int = 0, back_prop_type: int = 1, cropgrad=True, gradfreqs: numpy.ndarray = array([], shape=(1, 0), dtype=float64), param_type: int = 0, tmax: float = 0, tmin: float = 0, fmin: float = 0, fmax: float = 0, filter_offset: bool = False, offmin: float = -inf, offmax: float = inf, inputres: int = 0, restype: int = 0, scalerms: int = 0, scalermsnorm: int = 0, scaleshot: int = 0, seisout: int = 2, resout: int = 0, rmsout: int = 0, movout: int = 0, file: str = 'SeisCL', workdir: str = './seiscl') method of SeisCL.SeisCL.SeisCL instance
You can define parameters controlling forward and adjoint computations
in the __init__ method or afterward.
Parameters defining the spatial and temporal grid:
:param N: Grid size [NZ, NX] in 2D or [NZ, NY, NX] in 3D
:param ND: Flag for dimension.
3: 3D, 2: 2D P-SV, 21: 2D SH, 22: 2D acoustic
:param dh: Grid spatial spacing
:param dt: Time step size
:param NT: Number of time steps
Parameters for the Generalized Standard linear solid implementing
seismic attenuation.
:param L: Number of attenuation mechanism (L=0 elastic)
:param f0: Central frequency of the relaxation. Also used
as the peak frequency of the default ricker
wavelet source
:param FL: Frequencies of the attenuation mechanism
Parameters of the finite difference stencils
:param FDORDER: Order of the finite difference stencil
Values: 2, 4, 6, 8, 10, 12
:param MAXRELERROR: Select method to compute FD coefficients
0: Taylor-coeff.
Holberg-coeff with maximum phase velocity error
1: 0.1 % 2: 0.5 % 3: 1.0 % 4: 3.0 %
:param FP16: FP16 computation:
0:FP32
1:FP32 vectorized, (faster than 0)
2: FP16 IO (computation are performed in FP32)
3: FP16 IO + COMP (computation in FP16,
requires CUDA)
Parameters controlling the acquisition
:param src_pos_all: Position of all shots containted in dataset
Array [sx sy sz srcid src_type] x nb sources
The same srcid are fired simulatneously
src_type: 100: Explosive, 0: Force in X, 1:
Force in Y, 2:Force in Z
:param rec_pos_all: Position of all the receivers in the dataset
Array [gx gy gz srcid recid - - -] x nb traces
srcid is the source number
recid is the trace number in the record
:param src_all: Source signals. NT x number of sources
Parameters defining the Boundary conditions
:param freesurf: Include a free surface 0: no, 1: yes
:param abs_type: Absorbing boundary type:
1: CPML, 2: Absorbing layer of Cerjan
:param VPPML: Vp velocity near CPML boundary
:param NPOWER: Exponent used in CMPL frame update
:param FPML: Dominant frequency of the wavefield
:param K_MAX_CPML: Coefficient involved in CPML
(may influence simulation stability)
:param nab: Width in grid points of the absorbing layer
:param abpc: Exponential decay for absorbing layer of Cerjan
Parameters controlling parallelization
:param with_docker: If True, use a docker implementation of SeisCL
:param with_mpi: Use mpi parallelization (True) or not.
:param NP: Number of MPI processes to launch
:param pref_device_type: Type of processor used (OpenCL version):
2: CPU, 4: GPU, 8: Accelerator
:param MPI_NPROC_SHOT: Maximum number of MPI process (nodes or gpus)
involved in domain decomposition
:param nmax_dev: Maximum number of GPUs per process
:param no_use_GPUs: Array of device numbers that should not be used
Example: 1--On a node with 4 GPUs, to use all 4 GPUs, each gpu computes
a different shot: with_mpi=True, NP=4, nmax_dev=1
2-- On a node with 4 GPUS, all of them used for model
decomposition: with_mpi=True, NP=1, nmax_dev=4
Parameters controlling how the gradient is computed.
:param gradout: Output gradient 1:yes, 0: no
:param Hout: Output approximate Hessian 1:yes, 0: no
:param gradsrcout: Output source gradient 1:yes, 0: no
:param back_prop_type: Type of gradient calculation:
1: backpropagation (elastic only)
2: Discrete Fourier transform
:param cropgrad: If true, when calling read_grad, the gradient
in the absorbing boundaries are set to 0.
:param gradfreqs: Frequencies of gradient with DFT
:param param_type: Type of parametrization:
0:(rho,vp,vs,taup,taus)
1:(rho, M, mu, taup, taus),
2:(rho, Ip, Is, taup, taus)
:param tmax: Maximum time of gradient computation
:param tmin: Minimum time of gradient computation
:param fmin: Maximum frequency of gradient. A butterworh
filter is used to remove frequencies lower
than fmin
:param fmax: Minimum frequency of gradient. A butterworh
filter is used to remove frequencies lower
than fmin
:param filter_offset: If true, will only allow the minimum and
maximum offset during computation
:param offmin: Maximum offset to compute
:param offmax: Minimum offset to compute
:param inputres: If 1, the gradient computation needs the
residuals (adjoint sources) to be provided.
This allows to compute the cost and adjoint
sources in python, providing more flexibility.
If 0, SeisCL computes the cost and adjoint
sources.
The rest of these parameters apply if inputres=0
:param restype: Type of costfunction
0: l2 cost. 1: Cross-correlation of traces
:param scalerms: Scale each modeled and recorded traces
according to its rms value, then scale residual
by recorded trace rms when computing cost
:param scalermsnorm: Scale each modeled and recorded traces
according to its rms value before computing
cost
:param scaleshot: Scale all of the traces in each shot by the
shot total rms value when computing cost
Parameters setting the type of outputs
:param seisout: Output seismograms
1: output velocities,
2: output pressure,
3: output stresses, output everything
:param resout: Output residuals 1:yes, 0: no
:param rmsout: Output rms value of the cost 1:yes, 0: no
:param movout: Output movie every n frames
Parameters for file creation
:param file: Base name of the file to create for SeisCL.
Created files will be appended thr right suffix,
i.e. the model file will be file_model.mat
:param workdir: The name of the directory in which to create
the file for SeisCL
Simulation Grid
We begin setting the simulation grid, a 2D square domain of size \(\small 500 \times 500 \; m\) with a cell size of \(2m\) in all directions. We perform 1000 timesteps with a step size of with \(0.25ms\) for a simulation time of \(0.25s\).
[4]:
seis.ND = 2 # Number of dimension
seis.N = np.array([250, 250]) # Grid size [NZ, NX, NY]
seis.dh = dh = 2 # Grid spatial spacing
seis.dt = dt = 0.25e-03 # Time step size
seis.NT = NT = 1000 # Number of time steps
Velocity Model
For the elastic wave equation, three material parameters are required : \(V_p\), \(V_s\) and \(\rho\). For the sake of simplicity, we the parameters are homegeneous throughout the model, with \(V_p = 3500 m/s\), \(V_s = 2000 m/s\), \(\rho = 2000 kg/m^3\). We store all parameters in a dictionnary, which will be used by the SeisCL class.
[5]:
vp = 3500
vs = 2000
rho = 2000
vp_a = np.zeros(seis.csts['N']) + vp
vs_a = np.zeros(seis.csts['N']) + vs
rho_a = np.zeros(seis.csts['N']) + rho
model_dict = {"vp": vp_a, "rho": rho_a, "vs": vs_a}
Sources and receivers
Source and receivers are totally custumizable by SeisCL: each source can have different receivers, and multiple sources can be fired simultaneously. This flexibilty comes with a price: the definition of the source and receivers becomes verbose pretty easily.
We provide some utility methods to ease the definition of source and receivers. Here, we use the method SeisCL.surface_acquisition_2d, which create a surface acquisition geometry, with receivers and sources at a defined interval of the grid steps along all the model.
[6]:
seis.surface_acquisition_2d(ds=1000)
All this function did is to define the attributes SeisCL.src_pos_all and SeisCL.rec_pos_all.
[7]:
print(seis.src_pos_all)
[[ 36.]
[ 0.]
[ 36.]
[ 0.]
[ 100.]]
We see here that we have one source, at x=36m, z= 36m (just outside the absorbing boundary), and that the source type is 100 (explosive source). The format of this array is [x, y, z, srcid, src_type], where all shots having the same src_id are fired simultaneously.
[8]:
print(seis.rec_pos_all[:,0])
[ 40. 0. 36. 0. 1. 0. 0. 0.]
We see here that the first receiver is at x=40m and z=36m. The receiver is for the source id 0, and have the receiver id 1. The receiver id should be unique and start at 1 for the first receiver. The format of this array is [x, y, z, srcid, recid, src_type, none, none, none].
Simulation
we can now perform the simulation. The method SeisCL.set_forward is used to initialize the simulation, i.e. writing the hdf5 files containing simulation parameters and the model. The method SeisCL.execute() launches SeisCL (the compiled program) and wait for the end of execution.
[12]:
seis.set_forward([0], model_dict, withgrad=False)
seis.execute()
datafd = seis.read_data()[0]
All hdf5 files are created in a temporary directory, called seiscl by default
[13]:
os.listdir("./seiscl")
[13]:
['SeisCL_csts.mat', 'SeisCL_model.mat', 'SeisCL_cache', 'SeisCL_dout.mat']
Each file contains specific information required for simulations: * SeisCL_csts.mat contains all simulation parameters (dt, dh, …) * SeisCL_model.mat contains the material parameter arrays * SeisCL_dout.mat contains the simulated data.
Visualization
With matplotlib, we can finally analyze and visualize the result of our simulation.
[14]:
clip = 0.1
extent = [min(seis.rec_pos_all[0]), max(seis.rec_pos_all[0]), (datafd.shape[0]-1)*dt, 0]
vmax = np.max(datafd) * clip
vmin = -vmax
fig, ax = plt.subplots(1, 1, figsize=[4, 6])
ax.imshow(datafd, aspect='auto', vmax=vmax, vmin=vmin, extent = extent,
interpolation='bilinear', cmap=plt.get_cmap('Greys'))
ax.set_title("FD solution to elastic wave equation in 2D \n", fontsize=16, fontweight='bold')
ax.set_xlabel("Receiver position (m)")
ax.set_ylabel("time (s)")
plt.show()
Learn more
SeisCL has many more capabilities, including viscoelastic propagation in 2D and 3D, parallelization on multiples devices, adjoint modeling and reverse time migration. You can learn more on these different features with the following notebooks.
[ ]: