Model Interpolation

# initialize notebook
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import obspy
import salvus.namespace as sn

plt.rcParams["figure.figsize"] = (10, 8)

def get_marmousi():
    nx, ny = 2301, 751
    dx, dy = 4.0, 4.0

    rho = np.empty((ny, nx))
    st = obspy.read("data/marmousi_density.segy.gz")
    for _i, tr in enumerate(st):
        rho[_i, :] = tr.data

    vp = np.empty((ny, nx))
    st = obspy.read("data/marmousi_velocity.segy.gz")
    for _i, tr in enumerate(st):
        vp[_i, :] = tr.data

    x = np.arange(nx) * dx
    y = np.arange(ny) * dy

    return x, y, vp, rho

x, y, vp, rho = get_marmousi()
fig, axes = plt.subplots(2, 1)
axes[0].imshow(vp)
axes[1].imshow(rho)
plt.show()

fmax = 10.0  # maximum frequency in Hz
elements_per_wavelength = 2.0  # resolution criterion

# region in the marmousi model to mesh
mesh_min_x = x.max() * 0.2
mesh_max_x = x.max() * 0.7
mesh_min_y = y.min()
mesh_max_y = y.max()

diff_x = mesh_max_x - mesh_min_x
diff_y = mesh_max_y - mesh_min_y

m = sn.simple_mesh.basic_mesh.CartesianHomogeneousAcoustic2D(
    vp=vp.min(),
    rho=rho.min(),
    x_max=float(diff_x),
    y_max=float(diff_y),
    elements_per_wavelength=elements_per_wavelength,
    max_frequency=fmax,
    tensor_order=4,
).create_mesh()
m.points[:, 0] += mesh_min_x
m.points[:, 1] += mesh_min_y

from scipy.interpolate import RectBivariateSpline

f_vp = RectBivariateSpline(x, y, vp.T[:, ::-1])
f_rho = RectBivariateSpline(x, y, rho.T[:, ::-1])

nodes = m.get_element_nodes()

nodes_x = nodes[:, :, 0]
nodes_y = nodes[:, :, 1]

vp_nodes = f_vp(nodes_x, nodes_y, grid=False)
rho_nodes = f_rho(nodes_x, nodes_y, grid=False)

m.attach_field("VP", vp_nodes)
m.attach_field("RHO", rho_nodes)
m.attach_field("fluid", np.ones(m.nelem))
m