MondaicA common use case for spectral-element simulations is for the simulation of
wavefields in the presence of surface, subsurface, or seabed topography. In this tutorial we'll
- Define a set set of interfaces with topography by using interpolating splines.
- Generate a mesh which respects the defined topography.
- Set up and run an acoustic simulation through the resulting mesh.
Let's get started
by importing the relevant python dependencies.
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")The first thing we'll do is set up the splines which will define the
interfaces in our model. Before we get started though, let's set the x extent of our model, which will be km.
x_min = 0.0
x_max = 5000.0To help topography, we'll make use of the
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,
]We'll also define the p-wave velocities for each layer, and use some helper functions from Salvus Mesh to convert these values to realistic density and s-wave velocities.
# 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)The
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",
]Now we're ready to go ahead and generate our interpolating splines. Given the variables we defined above, we can call the function below to get the horizons of our model.
splines = sn.toolbox.get_interpolating_splines(
layers_x, layers_y, kind=interpolation_styles
)Below we simply just plot the splines to get an idea of what our model will look like in the end. Indeed the extents and horizons look as we expect, and we can move on to generating the mesh itself. Note that none of the horizons here cross each other, and if this were the case than the subsequent mesh generation would fail with an error. We'll re-visit how to deal with pinched-out layers and crossing horizons in a future tutorial.
# 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')
With our interfaces defined, we can now go ahead and generate the mesh itself. For this case a wrapper is included in the SalvusToolbox. This wrapper takes a series of parameters and attempts to build an optimal mesh with topography, and will also mark each individual region with a flag so it can be identified later. A more detailed description of the arguments one can pass is given below:
- x_min, x_max: The minimum and maximum x-values defining the model's horizontal extent. These are the dimensions to be considered before the attachment of absorbing boundaries. If absorbing boundaries are attached the resulting extent of the mesh will grow beyond x_min and x_max, but the solution will not be physical in the extruded regions.
- splines: The list of interpolation objects which we have generated above. These could of course be generated independently as well as long as their form is the same.
- absorbing_boundaries: A tuple consisting of the side sets which will be made absorbing, along with the number of elements to be extruded in each dimension for the absorbing layer. As previously stated, we recommend an absorbing layer of 3.5+ wavelengths if minimal reflections are desired.
- slowest_velocities: A list of the slowest velocities in each layer, ordered the same as the interfaces passed through splines (i.e. top to bottom). These values will be used, along with maximum_frequency and elements_per_wavelength to determine the element sizes in each region.
- use_refinements: Create a totally unstrucuted mesh by allowing for vertical refinements. Whether or not this will result in a performance increase is dependent on the model, we recommend try with and without refinements to see what works best for you.
- maximum_frequency: The maximum frequency to resolve in the mesh, which will be used in conjunction with elements_per_wavelength.
# 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),
)This function has returned a tuple of two elements. The first entry contains a list of mesh chunks, each of which correspond to a separate region in between two of the defined horizons. The second value returns the minimum size in meters of the extruded boundaries. This value will be used later on when running the actual simulation.
Our four individual mesh chunks are now properly sized and deformed, and are ready to be glued together. Each individual chunk has an elemental variable called
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)...and now the meshing stage is complete. The cell below will visualize the complete mesh. Note that we have the usual side sets (x0, x1, y0, y1) on the edges of our mesh, but we also have additional side sets which mark the internal discontinuities. These internal side sets can be handy when placing sources and receivers, as we'll explore in a future tutorial.
mesh # Visualize the mesh.
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7d8f90524090>
Attaching the parameters we defined above is simple because the elements are flagged by region.
nodes = mesh.get_element_nodes()[:, :, 0]
vp_a, vs_a, ro_a = np.ones((3, *nodes.shape))
for _i, (vp_val, vs_val, ro_val) in enumerate(zip(vp, vs, rho)):
# Find which elements are in a given region.
idx = np.where(mesh.elemental_fields["region"] == _i)
# Set parameters in that region to a constant value.
vp_a[idx] = vp_val
vs_a[idx] = vs_val
ro_a[idx] = ro_val
# Attach parameters.
for k, v in zip(["VP", "VS", "RHO"], [vp_a, vs_a, ro_a]):
mesh.attach_field(k, v)
# Attach acoustic / elastic flag.
mesh = sn.toolbox.detect_fluid(mesh)And of course we can visualize the mesh again to ensure everything is as we expect.
mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7d8f90524090>
This section has been convered in other tutorials, so we'll only briefly outline the main steps here.
First, we create a project based on the simulation domain defined above. We'll then add an event with some receivers placed directly on the surface. Note that because the surface is curved this process in non-trivial, and to assist we'll use the
SideSetVectorPoint2D class. Finally, we run a quick simulation and visualize the result in a shotgather.p = sn.Project.from_domain(
path="project",
domain=sn.domain.dim2.BoxDomain(x0=x_min, x1=x_max, y0=0.0, y1=2700.0),
load_if_exists=True,
)
p.viz.nb.domain()
Now, we add an event with a single source at the surface and an array of 1000 receivers 1m below that same deformed surface. At each receiver we record velocity and strain.
# Location of source (this will snap to the closet side-set).
loc = [x_max / 2, 2500.0]
# Create the source.
source = sn.simple_config.source.cartesian.SideSetVectorPoint2D(
fx=1,
fy=1,
point=loc,
direction="y",
side_set_name="y1",
)receivers = (
sn.simple_config.receiver.cartesian.SideSetHorizontalPointCollection2D(
x=np.linspace(x_min, x_max, 1000),
offset=-1.0,
side_set_name="y1",
station_code="xx",
fields=["velocity", "strain"],
)
)
p.add_to_project(
sn.Event(event_name="event", sources=source, receivers=receivers)
)Because we have created the mesh externally, it is added as a
UnstructuredMeshSimulationConfiguration to the project. Note that in this case it is the user's responsibility to ensure the mesh is accurate up to the desired frequency content of the source, and that boundary conditions have been properly set already.The
EventConfiguration is the same as for the ordinary SimulationConfiguration. After adding the UnstructuredMeshSimulationConfiguration to the project, it works like a simulation configuration, and we can refer to its name to launch simulations or query data.if "sim" not in p.simulations.list():
p.add_to_project(
sn.UnstructuredMeshSimulationConfiguration(
name="sim",
unstructured_mesh=mesh,
event_configuration=sn.EventConfiguration(
wavelet=sn.simple_config.stf.Ricker(
center_frequency=max_frequency / 2
),
waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
end_time_in_seconds=2.0
),
),
),
)When visualizing the simulation setup, we notice that the line of receivers is nicely following the surface topography.
p.viz.nb.simulation_setup("sim", events=["event"])
<salvus.flow.simple_config.simulation.waveform.Waveform object at 0x7d8f842e9790>
Now we are all set to call
launch and submit the simulation.p.simulations.launch(
simulation_configuration="sim",
events="event",
site_name=SALVUS_FLOW_SITE_NAME,
ranks_per_job=2,
)
p.simulations.query(block=True)[2024-11-20 13:06:58,622] INFO: Submitting job ... Uploading 1 files... 🚀 Submitted job_2411201306734055_6f6ca7a49c@local
True
Once the simulation is done (2D elastic, high accuracy up to 20 Hz -- should take less than 1 minute on 2 cores), we can then plot shotgathers of all the fields we've output (vector velocity, as well as the elastic strain).
waveform_data = p.waveforms.get(data_name="sim", events="event")[0]Let's start with velocity with components
X and Y.waveform_data.plot(
component="X",
receiver_field="velocity",
plot_types=["shotgather"],
)
waveform_data.plot(
component="Y",
receiver_field="velocity",
plot_types=["shotgather"],
)
We can also look at strain gathers. The components are enumerated (
0, 1, 2) in Voigt order.waveform_data.plot(
component="0",
receiver_field="strain",
plot_types=["shotgather"],
)
waveform_data.plot(
component="1",
receiver_field="strain",
plot_types=["shotgather"],
)
waveform_data.plot(
component="2",
receiver_field="strain",
plot_types=["shotgather"],
)
PAGE CONTENTS