A visual overview of the propagators

We provide in this notebook a comparison of the wavefields propagated with the different propagators implemented in th recipes. We consider a fixed anisotropic elastic multilayer model for all experiemnts.

import numpy as np
import matplotlib.pyplot as plt
import psutil

nthreads = int(np.ceil(psutil.cpu_count() / 3  * 2))

Utility functions

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

def plot_model(model):
    iy = param.shape[1]//2
    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[:, iy, :].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

# Available propagators
print("\n".join(recipes_registry.keys()))
Submodule devito version mismatch. Installed: 4.8.14+89.gd6980b7e9, Required: 4.8.14+86.g6d494dbdf

iso-acoustic
tti-zhang
tti-fletcher
tti-selfadjoint
iso-elastic
tti-elastic
vti-zhang
vti-fletcher
vti-selfadjoint
vti-elastic
# Setup grid
n = (64, 64, 64)
so = 8

args = {'time_order': 2, 'space_order': so, 'shape': n, 'fs': True,
        'nt' : 600, 'fd_kernel_only': False, 'density': True, 'abox': False,
        'dtype': np.float32, 'save': False, 'ndim': 3, 'nbl': 40, 'f0': 0.0075}

modelelas = layered_model('elastic-tti', nlayers=4, spacing=(20, 20, 20), **args)
modelacou = layered_model('acoustic-tti', nlayers=4, spacing=(20, 20, 20), dt=modelelas.critical_dt, **args)
modelelas.physical_params()
{'delta': delta(x, y, z),
 'eta': eta(x, y, z),
 'f': f(x, y, z),
 'mu': mu(x, y, z),
 'vp': vp(x, y, z),
 'theta': theta(x, y, z),
 'b': b(x, y, z),
 'epsilon': epsilon(x, y, z),
 'lam': lam(x, y, z),
 'phi': phi(x, y, z)}
plt.figure(figsize=(10, 10))
plt.subplot(4,2,1)
plt.imshow(modelelas.vp.data[:, 11, :].T, aspect="auto")
plt.title("P-wave velocity")
plt.colorbar()
plt.subplot(4,2,2)
plt.imshow(modelelas.mu.data[:, 11, :].T, aspect="auto")
plt.title(r"S-wave velocity (Lame \mu)")
plt.colorbar()
plt.subplot(4,2,3)
plt.imshow(modelelas.epsilon.data[:, 11, :].T, aspect="auto")
plt.title(r"$\epsilon$")
plt.colorbar()
plt.subplot(4,2,4)
plt.imshow(modelelas.delta.data[:, 11, :].T, aspect="auto")
plt.title(r"$\delta$")
plt.colorbar()
plt.subplot(4,2,5)
plt.imshow(modelelas.theta.data[:, 11, :].T, aspect="auto")
plt.title(r"$\theta$")
plt.colorbar()
plt.subplot(4,2,6)
plt.imshow(modelelas.phi.data[:, 11, :].T, aspect="auto")
plt.title(r"$\phi$")
plt.colorbar()
plt.tight_layout()

props = []
for k, prop in recipes_registry.items():
    if 'vti' in k:
        continue
    if 'elas' in k:
        model = modelelas
    else:
        model = modelacou
    props.append(prop(model, args))
ElasticTTI factorization: on-the-fly C_{ij}
props = sorted(props, key=lambda x: x.__class__.__name__)
from IPython.display import display, clear_output

plt.figure(figsize=(10, 10))

for (i, solver) in enumerate(props):
    name = solver.__class__.__name__
    src = solver.src
    src.coordinates.data[0, 0] /= 3
    if solver.time_order == 1:
        src.data[:, 0] = np.cumsum(props[0].src.data, axis=0)[:, 0] / props[0].model.critical_dt
    # get pressure
    solver.forward(save=False, nthreads=nthreads, src=src)
    p = np.array(solver.pressure_data)
    iy = p.shape[1] // 2
    c = np.quantile(np.abs(p[:, iy, :].reshape(-1)), .99)
    plt.subplot(3, 2, i+1)
    plt.imshow(p[:, iy, :].T, cmap="seismic", vmin=-c, vmax=c, aspect="auto")
    plt.title(name)

    display(plt.gcf())  # Show current figure
    clear_output(wait=True)  # Clear previous output (optional to avoid clutter)

plt.tight_layout()
# display(plt.gcf())  # Final display

Back to top