Mondaic - Full Waveform Solutions

This tutorial is presented as Python code running inside a Jupyter Notebook, the recommended way to use Salvus. To run it yourself you can copy/type each individual cell or directly download the full notebook, including all required files.

Download Jupyter Notebook [126.96 kB]

Copy

# Magic
%matplotlib inline
%config Completer.use_jedi = False

A 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.

Copy

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

5

km.

Copy

x_min = 0.0
x_max = 5000.0

To 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

x

and

y

. In the cell below we define a series of

5

discontinuities by specifying their interpolation points in each coordiante. We multiply by

1000

here simply for brevity.

Copy

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.

Copy

# 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.

Copy

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.

Copy

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.

Copy

# 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.

Copy

# 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. We can plot the mesh chunks individually using the build in plot functionality of the mesher. Keep in mind that this may take some time for bigger meshes, so if you are running at higher frequencies you may want to skip of comment out the following cell.

Copy

f, axs = plt.subplots(2, 2, figsize=(15, 10), sharex=True)
for _i, (ax, sub_mesh) in enumerate(zip(axs.flatten(), mesh)):
    plt.sca(ax)
    ax.set_xlabel("x (m)")
    ax.set_ylabel("y (m)")
    ax.set_title(f"Mesh Chunk {_i + 1}")
    sub_mesh.plot(show=False, figure=f, linewidths=0.1)

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...

Copy

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.

Copy

mesh  # Visualize the mesh.

<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fbaef653910>

Attaching the parameters we defined above is simple because the elements are flagged by region.

Copy

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.

Copy

mesh

<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fbaef653910>

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.

Copy

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.

Copy

# 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",
)

Copy

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.

Copy

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.

Copy

p.viz.nb.simulation_setup("sim", events=["event"])

<salvus.flow.simple_config.simulation.waveform.Waveform object at 0x7fbaef062c90>

Now we are all set to call launch and submit the simulation.

Copy

p.simulations.launch(
    simulation_configuration="sim",
    events="event",
    site_name=SALVUS_FLOW_SITE_NAME,
    ranks_per_job=2,
)
p.simulations.query(block=True)

[2024-03-15 09:07:02,237] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2403150907419315_9d712a347d@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).

Copy

waveform_data = p.waveforms.get(data_name="sim", events="event")[0]

Let's start with velocity with components X and Y.

Copy

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.

Copy

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"],
)