Mondaic
PricingAbout UsDocumentation
Mondaic
Version:
No results found.

Creating and plotting custom gathers from Salvus' receiver outputs

Salvus records the wavefields on receiver locations as specified by the user. These wavefields can be readily obtained for every receiver, but Salvus also has built-in functionality to retrieve data as well as plot sorted, filtered and selected subsets of these receivers.
This notebook illustrated how to use custom callback functions to create subsets of the data, and plot them.
Copy
import matplotlib.pyplot as plt
import numpy as np
import os
import xarray as xr
import dataclasses

import salvus.namespace as sn

# General configuration variables
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
PROJECT_DIR = "project"

Creating two synthetic datasets

To show how different receiver gathers can be plotted, we will need some data. In the following few cells we will create two synthetic datasets on a simple box domain, with very dense receiver spacing. The goal of this notebook is to show how to meaningfully collect data from such an array of receivers.

Creating a domain

We start with creating a box domain that is rather oblong.
domain = sn.domain.dim3.BoxDomain(
    x0=0, x1=0.1, y0=0.0, y1=0.1, z0=-0.0, z1=0.025
)

p = sn.Project.from_domain(PROJECT_DIR, domain, load_if_exists=True)

p.viz.nb.domain()
In this domain, let's add a source, and an array of receivers. The receivers aren't added exactly to the surface, but projected onto it from the top down, using SideSetPoint3D and with initial point [x, y, 2.0] and direction z, intersecting with side_set z1.
src_3d = sn.simple_config.source.cartesian.VectorPoint3D(
    x=0.03, y=0.065, z=0.025, fx=0.0, fy=0.0, fz=-1.0
)

# Create a receiver grid with regular spacing
x_locations_recs = [0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09]
y_locations_recs = np.linspace(0.01, 0.09, 50)

recs = [
    sn.simple_config.receiver.cartesian.SideSetPoint3D(
        point=[x, y, 2.0],
        direction="z",
        side_set_name="z1",
        station_code=f"TOP-X{int(ix)}-Y{int(iy)}",
        fields=["displacement"],
    )
    for ix, x in enumerate(x_locations_recs)
    for iy, y in enumerate(y_locations_recs)
]

# Package up sources and receivers into an event collection
p.add_to_project(
    sn.EventCollection.from_sources(
        sources=[src_3d],
        receivers=recs,
    )
)
Now, finally, let's add a model, waveform details and a simulation configuration, i.e., everything that is needed to create an appropriate mesh and run a simulation. We use a homogeneous background model, a fixed end time, and a delta function as a source time function.
# Create a homogeneous, isotropic, elastic model
bm = sn.model.background.homogeneous.IsotropicElastic(
    vp=1000.0, rho=1000.0, vs=1000 / (2.0**0.5)
)
mc_background = sn.ModelConfiguration(background_model=bm)

# Set a fixed end time, as predicted by the user
wsc = sn.WaveformSimulationConfiguration(end_time_in_seconds=0.00030)

# Use a maximum frequency of 50 kHz. This strongly influences the meshing
# resolution, and thereby the cost of running a simulation
max_frequency = 30e3

# Create an event configuration by setting the source time function
event_configuration = sn.EventConfiguration(
    waveform_simulation_configuration=wsc,
    wavelet=sn.simple_config.stf.Ricker(center_frequency=0.5 * max_frequency),
)

# Add everything into a simulation configuration and add it to the project
p.add_to_project(
    sn.SimulationConfiguration(
        name="simulation_homogeneous",
        elements_per_wavelength=2.0,
        tensor_order=1,
        max_frequency_in_hertz=max_frequency,
        model_configuration=mc_background,
        event_configuration=event_configuration,
    )
)
This is all we need to create synthetic data. However, it would be interesting later on to compare at least two different datasets. To this extent, let's create another simulation with a different model. This time, let's add some heterogeneity on top of our background model.
The following function creates an xarray that can be passed to sn.model.volume.cartesian.GenericModel.
def get_ellipsoidal_inclusion(
    domain,
    bm,
    vp_perturbation: float = 0.3,
    vs_perturbation: float = -0.2,
    rho_perturbation: float = 0.8,
):
    """
    Returns an xarray limited to the `domain` bounds. It centers a
    ellipsoidal inclusion of vp, vs and rho based on the values of the
    passed `bm`.

    """
    # Define resolution of xarray
    nx, ny, nz = 100, 100, 100

    # Get bounds
    (x0, x1), (y0, y1), (z0, z1) = domain.bounds.sorted_bounds

    # Get domain dimensions
    xl, yl, zl = x1 - x0, y1 - y0, z1 - z0

    # Get midpoint
    xm, ym, zm = (x1 + x0) / 2.0, (y1 + y0) / 2.0, (z1 + z0) / 2.0

    # Create coordinate arrays
    x = np.linspace(x0, y1, nx)
    y = np.linspace(y0, y1, ny)
    z = np.linspace(z0, z1, nz)
    xx, yy, zz = np.meshgrid(x, y, z, indexing="ij")

    # Create parameter background fields
    vp = bm.vp * np.ones_like(xx)
    vs = bm.vs * np.ones_like(xx)
    rho = bm.rho * np.ones_like(xx)

    # Compute ellipsoidal mask
    mask = (
        np.sqrt(
            ((xx - xm) / xl) ** 2
            + ((yy - ym) / yl) ** 2
            + ((zz - zm) / zl) ** 2
        )
        < 0.25
    )

    # Use the mask to perturb the model
    vp[mask] *= 1.0 + vp_perturbation
    vs[mask] *= 1.0 + vs_perturbation
    rho[mask] *= 1.0 + rho_perturbation

    # Convert the numpy arrays to an xarray dataset
    ds = xr.Dataset(
        data_vars={
            "vp": (["x", "y", "z"], vp),
            "vs": (["x", "y", "z"], vs),
            "rho": (["x", "y", "z"], rho),
        },
        coords={"x": x, "y": y, "z": z},
    )

    return ds


volumetric_model = get_ellipsoidal_inclusion(domain, bm)
Other than the fact that Salvus uses xarrays as a possible inputs for a volumetric model, they are also easily visualized. One can create a slice along a coordinate axis with .sel to verify our volume model.
plt.figure()
# Create a horizontal slice (z fixed to midpoint)
volumetric_model["vp"].sel(z=0.0125, method="nearest").plot()
plt.gca().set_aspect("equal")

plt.figure()
# Create a vertical slice (x fixed to midpoint). The transpose is taken simply
# to use the vertical axis where we expect it to.
volumetric_model["vs"].sel(x=0.05, method="nearest").transpose("z", "y").plot()
plt.gca().set_aspect("equal")

plt.figure()
# Create a vertical slice (x fixed to midpoint). The transpose is taken simply
# to use the vertical axis where we expect it to.
volumetric_model["rho"].sel(x=0.05, method="nearest").transpose(
    "z", "y"
).plot()
plt.gca().set_aspect("equal")
That model looks exactly like something interesting. Let's add it to the project, and create a new simulation the uses this model.
# Add the model with inclusion to the project
p.add_to_project(
    sn.model.volume.cartesian.GenericModel(
        name="model_with_inclusion", data=get_ellipsoidal_inclusion(domain, bm)
    )
)

# Create a model configuration to be used by the simulation configuration
mc_3d = sn.ModelConfiguration(
    background_model=None, volume_models="model_with_inclusion"
)

# Recreate the same simulation configuration, but with the new model
p.add_to_project(
    sn.SimulationConfiguration(
        name="simulation_with_inclusion",
        elements_per_wavelength=2.0,
        tensor_order=1,
        max_frequency_in_hertz=max_frequency,
        model_configuration=mc_3d,
        event_configuration=event_configuration,
    )
)
To verify our models, one might like to visualize the simulations. First the event without the inclusion.
p.viz.nb.simulation_setup("simulation_homogeneous", "event_0000")
[2024-07-05 10:09:38,740] INFO: Creating mesh. Hang on.