Freesurface for acoustic solvers

We provide in this notebook a comparison of the wavefields propagated with the different propagators implemented in the recipes that support a reflective free surface.

import numpy as np
import matplotlib.pyplot as plt

Utility functions

We setup a few utility function to plot the model and wavefield nicely

def plot_model(model):
    domain_size = 1.e-3 * np.array(param.grid.domain_size)
    extent = [model.origin[0], model.origin[0] + domain_size[0],
              model.origin[1] + domain_size[1], model.origin[1]]
    plt.imshow(param.data.T, aspect="auto", cmap="vik", extent=extent)
    plt.title("%s" % param.name)
    plt.colorbar()

Setup model

We first setup the model using the most complex physics available tti-elastic. This model contains all necessary parameters to propagate waves with all propagators.

from recipes import recipes_registry, layered_model
# Setup grid
n = (201, 201)
so = 8

args = {'time_order': 2, 'space_order': so, 'shape': n,
        'nt' : 250, 'fd_kernel_only': False, 'density': False, 'abox': False,
        'dtype': np.float32, 'save': False, 'ndim': 3, 'nbl': 20}

model = layered_model('elastic-tti', nlayers=4, fs=True, **args)
Operator `initdamp` ran in 0.01 s
model.physical_parameters
('damp', 'eta', 'vs', 'delta', 'epsilon', 'f', 'vp', 'b', 'theta')
plt.figure(figsize=(10, 10))
plt.subplot(4,2,1)
plt.imshow(model.vp.data.T, aspect="auto")
plt.title("P-wave velocity")
plt.colorbar()
plt.subplot(4,2,2)
plt.imshow(model.vs.data.T, aspect="auto")
plt.title("S-wave velocity")
plt.colorbar()
plt.subplot(4,2,3)
plt.imshow(model.epsilon.data.T, aspect="auto")
plt.title(r"$\epsilon$")
plt.colorbar()
plt.subplot(4,2,4)
plt.imshow(model.delta.data.T, aspect="auto")
plt.title(r"$\delta$")
plt.colorbar()
plt.subplot(4,2,5)
plt.imshow(model.theta.data.T, aspect="auto")
plt.title(r"$\theta$")
plt.colorbar()
plt.tight_layout()

props = [recipes_registry[k](model, args) for k in ['iso-acoustic', 'tti-zhang']]

Lets first show the clear reflection from the top

for solver in props:
    solver.src.coordinates.data[0, -1] = 200
    solver.forward(**args)
Operator `ForwardAcousticIsotropic` ran in 0.02 s
Operator `ForwardZhangTTI` ran in 0.05 s
plt.figure(figsize=(10, 10))

for (i, solver), name in zip(enumerate(props), recipes_registry):
    # get pressure
    p = np.array(solver.pressure_data)
    c = np.quantile(np.abs(p.reshape(-1)), .98)
    plt.subplot(1, 2, i+1)
    plt.imshow(p.T, cmap="seismic", vmin=-c, vmax=c, aspect="auto")
    plt.title(name)

plt.tight_layout()

And with a shallow source to see the phae effect

for solver in props:
    solver.src.coordinates.data[0, -1] = 10
    solver.forward(**args)
Operator `ForwardAcousticIsotropic` ran in 0.03 s
Operator `ForwardZhangTTI` ran in 0.05 s
plt.figure(figsize=(10, 10))

for (i, solver), name in zip(enumerate(props), recipes_registry):
    # get pressure
    p = np.array(solver.pressure_data)
    c = np.quantile(np.abs(p.reshape(-1)), .98)
    plt.subplot(1, 2, i+1)
    plt.imshow(p.T, cmap="seismic", vmin=-c, vmax=c, aspect="auto")
    plt.title(name)

plt.tight_layout()

Back to top