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 [127.31 kB]

The AGL Elastic Marmousi model is the most recent generation of the classic Marmousi seismic benchmark. In this tutorial we will cover how one can use this model within Salvus. In particular, we will cover:

How to read a SEG-Y file using Obspy
How to read in a rectilinear model using xarray
How to interpolate a 2-D regularly-gridded model onto the mesh
How to generate a synthetic seismic shotgather using the model

As always, our first job is to import the relevant Python packages into our notebook.

%matplotlib inline
%config Completer.use_jedi = False

import numpy as np
import os
import shutil
import pathlib
import requests
import time
import xarray as xr

import salvus.namespace as sn
from salvus.namespace import simple_config as sc


SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
SIMULATION_MAX_FREQUENCY = 5.0

Information about the Marmousi model we'll be using can be found here. The model is stored as a SEG-Y file, and the next cell will download and unpack the model for you. Note that the model file is relatively large (~150 MB), so make sure you have enough disk space before continuing.

# Download Marmousi model ~150 MB.
target_path = "elastic-marmousi-model.tar.gz"
if not pathlib.Path(target_path).exists():
    url = (
        "https://s3.amazonaws.com/open.source.geoscience/open_data"
        "/elastic-marmousi/elastic-marmousi-model.tar.gz"
    )

    response = requests.get(url, stream=True)
    if response.status_code == 200:
        with open(target_path, "wb") as f:
            f.write(response.raw.read())
        shutil.unpack_archive(target_path)

With our environment is set up, and the model downloaded, we can begin building the Salvus mesh. Most material models in seismic exploration are discretized onto a regular grid, and the open-source SalvusPythonToolkit provides a simple interface to create two or three dimensional Salvus meshes based on such models. As models can be stored in a wide variety of format such as SEG-Y, or RSF, we use xarray as an intermediary to encapsulate both the model parameters and geometry. More info on xarray, including extensive documentation, can be found here. Reading our Marmousi model into an xarray.Dataset is trivial, and the process can be inspected by opening the acompanying getting_started_tools.py file.

In general, if you are working working with a regularly gridded model in either two or three dimensions, a recipie similar to the function provided should be all you need to do to get your model into a form Salvus Mesh can understand. All that is required is an xarray dataset with the following information

The coordinate values in meters (named 'x', 'y', and optionally 'z')
The material parameters in SI units, with names following the parameter naming convention

# Read and plot marmousi model.
model_directory = pathlib.Path("./elastic-marmousi-model/model")
marmousi = sn.toolbox.read_elastic_marmousi(
    model_directory, ["VP", "RHO", "VS"]
)

A nice feature of xarray.Dataset object is the ability to quickly plot the contained data. Let's do this below as a way to quickly QC that the model reading went correctly.

marmousi["VP"].transpose("y", "x").plot(aspect="auto", figsize=(15, 5))
marmousi["VS"].transpose("y", "x").plot(aspect="auto", figsize=(15, 5))
marmousi["RHO"].transpose("y", "x").plot(aspect="auto", figsize=(15, 5))

<matplotlib.collections.QuadMesh at 0x7f13434dff10>

Note that the top layer models a fluid ocean. In order to accurately incorporate the interface of the ocean bottom into the spectral-element mesh. We will first cut out this part from the model and later on add it back as an ocean layer. This will ensure that the interface conditions between solid and fluid media are correctly imposed and that the model is indeed discontinuous across this interface.

marmousi_elastic = marmousi.where(marmousi["y"] <= 3500 - 37 * 12.5, drop=True)
marmousi_elastic["VP"].transpose("y", "x").plot(aspect="auto", figsize=(15, 5))
marmousi_elastic["VS"].transpose("y", "x").plot(aspect="auto", figsize=(15, 5))
marmousi_elastic["RHO"].transpose("y", "x").plot(
    aspect="auto", figsize=(15, 5)
)

<matplotlib.collections.QuadMesh at 0x7f1343a28210>

Yup, all looks good!

We'll use SalvusProject to manage our workflows. First, we initialize the project (or load it if it is already on disk). Since our simulation domain is defined by the extent of our Marmousi model, we'll use the from_volume_model constructor to initialize our project.

project_folder = pathlib.Path("project")

# Initialize a volumetric model from an xarray dataset.
mt = sn.model.volume.cartesian.GenericModel(
    name="marmousi_elastic", data=marmousi_elastic
)

# Create a project if it does not yet exist or load it from disk.
p = sn.Project.from_volume_model(
    path=project_folder, volume_model=mt, load_if_exists=True
)

Up to now, the domain only contains the solid parts, hence its vertical extent is only 3037.5 m.

p.viz.nb.domain()

Next, we also add a purely acoustic version of the Marmousi model to the project. Again, we explicitly model the upper part as a fluid ocean layer to make sure the velocity model is not smeared along the ocean bottom interface.

marmousi = sn.toolbox.read_elastic_marmousi(model_directory, ["VP", "RHO"])
marmousi_acoustic = marmousi.where(
    marmousi["y"] <= 3500 - 37 * 12.5, drop=True
)
p.add_to_project(
    sn.model.volume.cartesian.GenericModel(
        name="marmousi_acoustic", data=marmousi_acoustic
    )
)

Events in SalvusProject are containers to store waveform data. They are associated with sources and receivers.

Here we'll choose a scalar source type (as we're working in an acoustic medium), and use the side-set interface to allow for easy and intuitive source placement. Note that, since this is a purely acoustic simulation, if we chose an elastic source type (either Vector or MomentTensor) the simulation would give an error explaining that the given source type could not be attached to any element.

Note that at this stage, we only include the spatial characterization of the source. The source time function will enter later, once we specify the frequency of the simulation.

source = sn.simple_config.source.cartesian.ScalarPoint2D(
    x=8500.0, y=3490.0, f=1
)

Setting up receivers is similar to setting up sources. For this example we want to use a regularly-spaced grid of receivers which span the length of the domain. For this, we can use the SideSetHorizontalPointCollection2D receiver type. This allows us to choose a side-set (here, "y1"), and to place an array of receivers a certain distance from this side-set. In this case we choose to place the receivers righ at the ocean bottom. We also need to specify an array of dynamic fields to record. Here we choose pressure, which is equal to "phi_tt".

receivers = (
    sn.simple_config.receiver.cartesian.SideSetHorizontalPointCollection2D(
        x=np.linspace(0.0, 17000.0, 1000),
        offset=0.0,
        station_code="xx",
        side_set_name="y1",
        fields=["phi_tt"],
    )
)

Because both the source and the receivers are outside of the elastic domain, we skip the consistency check in the following, which would otherwise complain the vertical location.

p.events.add(
    sn.EventCollection.from_sources(sources=source, receivers=receivers),
    check_consistency_with_domain=False,
)

To set up our simulation, we need to specify a few more parameters. The first thing we'll do is set the duration of the simulation, by specifying end_time_in_seconds. We could additionally provide parameters like time_step_in_seconds and start_time_in_seconds. If they are not set, Salvus will choose sensible defaults for all three of these parameters. The start and end times will default to the values needed to smoothly simulate the entire source wavelet, and the time step will be computed inside Salvus as per the CFL criterion. Here we set the end time to 8 seconds, and keep the rest of the parameters at their default values.

wsc = sn.WaveformSimulationConfiguration(end_time_in_seconds=8.0)

We now have enough information to describe an EventConfiguration, which brings together the spatial and temporal aspects of our simulation. We will be using a Ricker wavelet for this simulation, with a center frequency of half the maximum frequency we've set for the mesh.

wavelet = sn.simple_config.stf.Ricker(
    center_frequency=SIMULATION_MAX_FREQUENCY / 2.0
)
ec = sn.EventConfiguration(
    waveform_simulation_configuration=wsc, wavelet=wavelet
)

Now we prepare for adding the ocean layer. The ocean bottom is flat, so the bathymetry model is just a straight line. The width of the ocean column corresponds exactly to the part we cut out at the beginnning of the tutorial.

x = np.linspace(-1e6, 1e6, 10)
y = 37 * 12.5 * np.ones_like(x)
ds = xr.Dataset(
    data_vars={"dem": (["x"], y, {"reference_elevation": 0.0})},
    coords={"x": x},
)
p.add_to_project(
    sn.bathymetry.cartesian.OceanLayer(
        name="ocean",
        data=ds,
        reference_elevation=0.0,
        ocean_surface_datum=3500.0,
    ),
    overwrite=True,
)

We'll also need to specify the absorbing boundary parameters we set above to mock an infinite domain extending from the other side-sets.

# Define coupled Clayton-Enqist / Kosloff
# damping boundaries at the the remaining side-sets.
abc = sn.AbsorbingBoundaryParameters(
    reference_velocity=2000.0,
    number_of_wavelengths=3.5,
    reference_frequency=wavelet.center_frequency,
)

With all of our components defined, we can now instantiate a concrete simulation configuration...

p.add_to_project(
    sn.SimulationConfiguration(
        name="marmousi_elastic",
        event_configuration=ec,
        absorbing_boundaries=abc,
        elements_per_wavelength=1.0,
        tensor_order=4,
        max_frequency_in_hertz=float(SIMULATION_MAX_FREQUENCY),
        model_configuration=sn.ModelConfiguration(
            background_model=None, volume_models="marmousi_elastic"
        ),
        bathymetry_configuration=sn.BathymetryConfiguration("ocean"),
    ),
    overwrite=False,
)

...and then view it below with all the events overlain. Hovering over each source or receiver will reveal more details.

p.viz.nb.simulation_setup("marmousi_elastic", events=p.events.list())

[2024-03-15 09:08:56,750] INFO: Creating mesh. Hang on.

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

p.add_to_project(
    sn.SimulationConfiguration(
        name="marmousi_acoustic",
        event_configuration=ec,
        absorbing_boundaries=abc,
        elements_per_wavelength=1.0,
        tensor_order=4,
        max_frequency_in_hertz=float(SIMULATION_MAX_FREQUENCY),
        model_configuration=sn.ModelConfiguration(
            background_model=None, volume_models="marmousi_acoustic"
        ),
        bathymetry_configuration=sn.BathymetryConfiguration("ocean"),
    ),
    overwrite=False,
)

p.viz.nb.simulation_setup("marmousi_acoustic", events=p.events.list())

[2024-03-15 09:09:02,134] INFO: Creating mesh. Hang on.

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

Next, we finally launch the simulation in the standard way as below. Remember: if the site you pass is remote, all the shuttling of data to/from the remote machine will be handled for you. Note: since we are using the SalvusProject framework, the simulations will only be submitted if they have not already been. This means that if you run the notebook multiple times, no redundant computations will take place.

The commented lines with the optional argument extra_output_configuration show how to specify volumetric (2- or 3-D) output for wavefield visualization and analysis. We don't output these by default because the file sizes involved can be quite large (in this case, the parameters below will produce 375 MB of output). If your machine has the capacity, feel free to comment the following in to generate time-varying volumetric output which can be visualized in Paraview.

for sim in ["marmousi_elastic", "marmousi_acoustic"]:
    p.simulations.launch(
        events=p.events.list(),
        site_name=SALVUS_FLOW_SITE_NAME,
        ranks_per_job=4,
        simulation_configuration=sim,
        # # uncomment for volume output
        # extra_output_configuration={
        #     "volume_data": {
        #         "sampling_interval_in_time_steps": 100,
        #         "fields": ["phi_tt", "stress"],
        #     },
        # }
    )

[2024-03-15 09:09:04,362] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2403150909544864_cd6f21968e@local

[2024-03-15 09:09:05,105] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2403150909136502_2cdce81e97@local

The previous cell launched our simulations and immediately returned, i.e. it did not wait for the simulations to complete. This is by design, and allows one the flexibility to shutdown and restart the notebook without loosing one's progress. In this case, however, we do want to wait until the simulations are complete before moving on.

To ensure that we block, we pass an optional arguemnt to query(). This will show update the progress info and wait until we can safely continue with the simulated data. If SalvusProject detects that the simulations have already been run and the results have been stored, the next cell is a no-op.

p.simulations.query(block=True)

True

Once the simulation has finished, the simulation results will automatically be added to the project.

We can query waveforms for each simulation configuration and event and directly visualize them in the notebook.

p.waveforms.get(data_name="marmousi_elastic", events=p.events.list())[0].plot(
    component="A",
    receiver_field="phi_tt",
    plot_types=["shotgather", "wiggles"],
)

p.waveforms.get(data_name="marmousi_acoustic", events=p.events.list())[0].plot(
    component="A",
    receiver_field="phi_tt",
    plot_types=["shotgather", "wiggles"],
)

Looks pretty cool!

Heterogeneous models and slice extraction

3D surface topography

© Mondaic AG (2024)

Site Map Contact Us

Academic License Agreement

PAGE CONTENTS