This documentation is not for the latest stable Salvus version.
# initialize notebook
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import obspy
from salvus.mesh import mesh_block
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()
hmax = vp.min() / fmax / elements_per_wavelength
nelem_x = int((mesh_max_x - mesh_min_x) / hmax) + 1
nelem_y = int((mesh_max_y - mesh_min_y) / hmax) + 1
sg = mesh_block.generators.cartesian.rectangle_2d(
nelem_x=nelem_x,
nelem_y=nelem_y,
min_x=mesh_min_x,
max_x=mesh_max_x,
min_y=mesh_min_y,
max_y=mesh_max_y,
)
m = sg.get_unstructured_mesh()
m.change_tensor_order(4)
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