Wave propogation in presence of obstacles

[nicl-nno]

Hello,

Thank you for the interesting project!

I have a question: is it possible to take into account the variable-shaped obstacles during the simulation of wave propagation with simwave?

E.g. as the figure from your paper presents:

?

[joao-bapdm]

Hi,

Thanks Nikolay,

If I understand the question correctly, arbitrary obstacles can be added to the model.

For instance, in Listing 1 the P-wave Marmousi velocity model available from https://s3.amazonaws.com/open.source.geoscience/open_data/elastic-marmousi/elastic-marmousi-model.tar.gz is used. The function read_2D_segy() loads the data as a numpy array, which can be manipulated as needed. If one would want to add a structure with anomalous high velocities to the model, one way to do it is:

from simwave import (
    SpaceModel, TimeModel, RickerWavelet, Solver, Compiler,
    Receiver, Source, read_2D_segy,
    plot_wavefield, plot_shotrecord, plot_velocity_model
)

import numpy as np
from skimage.draw import polygon


# Marmousi2 velocity model
marmousi_model = read_2D_segy('MODEL_P-WAVE_VELOCITY_1.25m.segy')

# Add obstacle
spacing = 30
orig_z, orig_x = 1930, 9870
poly = np.array((
    (orig_z + 0 * spacing, orig_x + 0 * spacing),
    (orig_z + 8 * spacing, orig_x + 4 * spacing),
    (orig_z + 12 * spacing, orig_x + 4 * spacing),
    (orig_z + 16 * spacing, orig_x + 8 * spacing),
    (orig_z + 20 * spacing, orig_x + 12 * spacing),
    (orig_z + 20 * spacing, orig_x + 16 * spacing),
    (orig_z + 16 * spacing, orig_x + 16 * spacing),
    (orig_z + 12 * spacing, orig_x + 20 * spacing),
    (orig_z + 16 * spacing, orig_x + 24 * spacing),
    (orig_z + 16 * spacing, orig_x + 28 * spacing),
    (orig_z + 16 * spacing, orig_x + 32 * spacing),
    (orig_z + 12 * spacing, orig_x + 28 * spacing),
    (orig_z + 8 * spacing, orig_x + 24 * spacing),
    (orig_z + 4 * spacing, orig_x + 20 * spacing),
    (orig_z + 4 * spacing, orig_x + 16 * spacing),
    (orig_z + 0 * spacing, orig_x + 12 * spacing),
    (orig_z + 0 * spacing, orig_x + 0 * spacing),
))
poly_rows, poly_columns = polygon(
    poly[:,0], poly[:,1], marmousi_model.shape
)
marmousi_model[poly_rows, poly_columns] = marmousi_model.max() + 800

compiler = Compiler(
    cc='gcc',
    language='cpu_openmp',
    cflags='-O3 -fPIC -ffast-math -std=c99'
)

# create the space model
space_model = SpaceModel(
    bounding_box=(0, 3500, 0, 17000),
    grid_spacing=(10.0, 10.0),
    velocity_model=marmousi_model,
    space_order=4,
    dtype=np.float64
)

# config boundary conditions
space_model.config_boundary(
    damping_length=(0, 700, 700, 700),
    boundary_condition=(
        "null_neumann", "null_dirichlet",
        "null_dirichlet", "null_dirichlet"
    ),
    damping_polynomial_degree=3,
    damping_alpha=0.001
)

# create the time model
time_model = TimeModel(
    space_model=space_model,
    tf=2.0,
    saving_stride=0
)

# create the set of sources
source = Source(
    space_model,
    coordinates=[(20, 8500)],
    window_radius=1
)

# crete the set of receivers
receiver = Receiver(
    space_model=space_model,
    coordinates=[(20, i) for i in range(0, 17000, 10)],
    window_radius=1
)

# create a ricker wavelet with 10hz of peak frequency
ricker = RickerWavelet(10.0, time_model)

# create the solver
solver = Solver(
    space_model=space_model,
    time_model=time_model,
    sources=source,
    receivers=receiver,
    wavelet=ricker,
    compiler=compiler
)

# run the forward
u_full, recv = solver.forward()

# remove damping extension from u_full
u_full = space_model.remove_nbl(u_full)

extent = [0, 17000, 3500, 0]

# plot the velocity model
plot_velocity_model(space_model.velocity_model, extent=extent)

# plot the last wavefield
plot_wavefield(u_full[-1], extent=extent)

# plot the seismogram
plot_shotrecord(recv)

Where the modified velocity model becomes:

For convenience I used the scikit-image package to create the obstacle above.

[nicl-nno]

Thank you for the example.
I'll try it.