import os
import matplotlib.pyplot as plt
import numpy as np
import salvus.namespace as sn
from salvus.mesh.simple_mesh import rho_from_gardners, vs_from_poisson
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
x_min = 0.0
x_max = 5000.0
get_interpolating_splines
function from SalvusToolbox
. This function
requires a set of interpolation points which can be used to anchor the
splines, and these points in turn are defined by a set of locations in
and . In the cell below we define a series of discontinuities by specifying
their interpolation points in each coordiante. We multiply by here simply for brevity.layers_x = [
np.array([0.0, 0.2, 1.0, 2.0, 3.0, 4.0, 4.8, 5.0]) * 1000,
np.array([0.0, 0.2, 1.0, 2.0, 3.0, 4.0, 4.8, 5.0]) * 1000,
np.array([0.0, 0.2, 1.0, 2.0, 3.0, 4.0, 4.8, 5.0]) * 1000,
np.array([0.0, 1.5, 3.5, 5.0]) * 1000,
np.array([0.0, 2.5, 5.0]) * 1000,
]
layers_y = [
np.array([2.0, 2.0, 1.9, 1.7, 2.0, 2.1, 2.0, 2.0]) * 1000,
np.array([1.6, 1.6, 1.5, 1.4, 1.3, 1.4, 1.5, 1.5]) * 1000,
np.array([0.5, 0.5, 0.7, 0.6, 1.1, 0.9, 1.2, 1.2]) * 1000,
np.array([0.2, 0.2, 0.4, 0.4]) * 1000,
np.array([0.0, 0.0, 0.0]) * 1000,
]
# Define p-velocities.
vp = np.array([2000.0, 2500.0, 2800.0, 3200.0])
# Compute vs and rho.
vs = vs_from_poisson(vp)
rho = rho_from_gardners(vp)
get_interpolating_splines
function also accepts an array of interpolations methods. Internally the function uses interp1d
from Scipy, so please check out the documentation here for a complete list of interpolation styles.interpolation_styles = [
"quadratic",
"quadratic",
"quadratic",
"linear",
"linear",
]
splines = sn.toolbox.get_interpolating_splines(
layers_x, layers_y, kind=interpolation_styles
)
# Plot the interfaces.
f = plt.figure(figsize=(10, 5))
x_plot = np.linspace(x_min, x_max)
for top, bot in splines:
plt.plot(x_plot, top(x_plot))
plt.plot(x_plot, bot(x_plot))
plt.xlabel("x (m)")
plt.ylabel("y (m)")
plt.title("Interfaces")
Text(0.5, 1.0, 'Interfaces')
# Maximum frequency to resolve with elements_per_wavelength.
max_frequency = 20.0
# Generate the mesh
mesh, bnd = sn.toolbox.generate_mesh_from_splines_2d(
x_min=0,
x_max=x_max,
splines=splines,
elements_per_wavelength=2,
maximum_frequency=max_frequency,
use_refinements=True,
slowest_velocities=vs,
absorbing_boundaries=(["x0", "x1", "y0"], 10.0),
)
region
attached to it, and we'll use this variable later on to track which region is which when we're attaching parameters. To sum the mesh regions into one continuous mesh all we need to do is simply call the sum
operation on the array of sub-meshes as below...mesh = np.sum(mesh)
# Add info about absorbing boundaries
mesh.attach_global_variable("max_dist_ABC", bnd)
mesh.attach_global_variable("ABC_side_sets", ", ".join(["x0", "x1", "y0"]))
mesh.attach_global_variable("ABC_vel", min(vs))
mesh.attach_global_variable("ABC_freq", max_frequency / 2.0)
mesh.attach_global_variable("ABC_nwave", 5.0)
mesh # Visualize the mesh.