Mondaic - Full Waveform Solutions

(Full)-Waveform Inversion (FWI) can reveal the properties of complex media that are otherwise not accessible to direct observation. This is based on measurements of mechanical waves - excited by either active or passive sources - that propagate through an object of interest and for which we record time series at certain receiver locations. Those data contain valuable information about the object's interior that we can use to create quantitative images of its material properties.

In this tutorial, we consider a fairly simple setup in 2D, which is inspired by typical apertures in ultrasound computed tomography (USCT) to illuminate human tissue.

Copy

%matplotlib inline

# This notebook will use this variable to determine which
# remote site to run on.
import os

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

Copy

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import salvus.namespace as sn

This is a purely synthetic study, so we have to generate the "measurements" ourselves using a model with a few inclusions. We use a simple box domain of 20 x 20 cm centered around the origin and insert an object with two spherical inclusions that mimics a breast phantom in a USCT experiment. The model for density (RHO) and speed of sound (VP) is created from a structured grid.

Copy

def get_spherical_inclusion():
    nx, ny = 200, 200
    x = np.linspace(-0.1, +0.1, nx)
    y = np.linspace(-0.1, +0.1, nx)
    xx, yy = np.meshgrid(x, y, indexing="ij")

    # Add 3 spherical inclusions
    vp = 1500.0 * np.ones_like(xx)
    rho = 990.0 * np.ones_like(xx)
    mask = np.sqrt(xx**2 + yy**2) < 0.05
    vp[mask] = 1480.0
    rho[mask] = 975.0

    mask = np.sqrt(xx**2 + (yy - 0.025) ** 2) < 0.015
    vp[mask] = 1555.0
    rho[mask] = 1040.0

    mask = np.sqrt(xx**2 + (yy + 0.025) ** 2) < 0.015
    vp[mask] = 1460.0
    rho[mask] = 1010.0

    ds = xr.Dataset(
        data_vars={
            "vp": (["x", "y"], vp),
            "rho": (["x", "y"], rho),
        },
        coords={"x": x, "y": y},
    )

    return ds

Let's take a look at the model.

Copy

true_model = get_spherical_inclusion()

# Plot the xarray dataset.
plt.figure(figsize=(16, 6))
plt.subplot(121)
true_model.vp.T.plot()
plt.subplot(122)
true_model.rho.T.plot()
plt.suptitle("Model with spherical inclusions")
plt.show()

This is also the starting point for setting up a project for the inversion. We can directly create the project from the model defined above. We will call this model the true_model.

Copy

# Uncomment the following line to delete a
# potentially existing project for a fresh start
# !rm -rf project
vm = sn.model.volume.cartesian.GenericModel(name="true_model", data=true_model)
p = sn.Project.from_volume_model(
    path="project",
    volume_model=vm,
    load_if_exists=True,
)

We assume a ring-shaped aperture with ultrasound transducers surrounding the breast phantom. To keep the computations cheap, we will use only 5 emitters and 100 receiving transducers which are the same for every emitter. The simple_config module has a few built-in options to create lists of sources and receivers, which we want to make use of. By defining the two rings below - one for sources and one for the receivers, we can readily create an EventCollection, i.e., a number of experiments characterized by the locations of the emitting and the receiving transducers.

Copy

srcs = sn.simple_config.source.cartesian.collections.ScalarPoint2DRing(
    x=0, y=0, radius=0.09, count=5, f=1.0
)

for _i, s in enumerate(srcs._sources):
    all_recs = sn.simple_config.receiver.cartesian.collections.RingPoint2D(
        x=0, y=0, radius=0.09, count=100, fields=["phi"]
    )
    recs = [
        r
        for r in all_recs._receivers
        if np.sqrt(
            (s.location[0] - r.location[0]) ** 2
            + (s.location[1] - r.location[1]) ** 2
        )
        > 0.03
    ]
    p.add_to_project(
        sn.EventCollection.from_sources(
            sources=[s], receivers=recs, event_name_starting_index=_i
        )
    )

The transducers are identical except for their location, so we use the same simulation time of 1.5 ms (and later on also the same source time function) for each of them. The WaveformSimulationConfiguration would also be the place to set the time step explicitly. Note, however, that Salvus determines the time step based on the Courant-Friedrich-Levy (CFL) condition if not provided explicitly. In the context of inversions, Salvus will automatically interpolate simulated and observed waveforms onto the same temporal grid before computing the misfit, so we do not need to worry about it. That's convenient, isn't it?

Copy

wsc = sn.WaveformSimulationConfiguration(end_time_in_seconds=0.00015)

Now, we put everything together to configure the simulations for the ground truth. To keep the costs low, we only consider a center frequency of 50 kHz and a mesh accurate up to 100 kHz.

Copy

ec = sn.EventConfiguration(
    waveform_simulation_configuration=wsc,
    wavelet=sn.simple_config.stf.Ricker(center_frequency=50000.0),
)
p.add_to_project(
    sn.SimulationConfiguration(
        name="true_model_100kHz",
        #
        # Settings that influence the mesh.
        elements_per_wavelength=2,
        tensor_order=4,
        max_frequency_in_hertz=100_000.0,
        #
        model_configuration=sn.ModelConfiguration(
            background_model=None, volume_models="true_model"
        ),
        # Potentially event dependent settings.
        event_configuration=ec,
    )
)

Now it is time for the first simulation in our project. This is only to create data though. We haven't even started with the inversion yet... Before running the simulation, let us have a look at the experimental setup.

Copy

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

[2025-07-23 21:52:54,976] INFO: Creating mesh. Hang on.

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

Now, we can launch the simulations.

Copy

p.simulations.launch(
    simulation_configuration="true_model_100kHz",
    events=p.events.list(),
    site_name=SALVUS_FLOW_SITE_NAME,
    ranks_per_job=1,
)

[2025-07-23 21:52:55,443] INFO: Submitting job array with 5 jobs ...

Copy

p.simulations.query(block=True)

True

Once the simulations are finished, we can query the data and visualize it directly in the notebook. Note that you might have to execute this cell again, in case the simulations were not already finished.

Copy

true_data = p.waveforms.get(
    data_name="true_model_100kHz",
    events=p.events.list(),
)

Similar to the event geometry, we can now slide through the individual recordings of each event.

Copy

true_data[0].plot(component="A", receiver_field="phi")

This concludes the first part of the tutorial. In the next step, we start with an initial guess of the model and compare the simulated data to the ones from the true model.

Academic License Agreement

Credits

PAGE CONTENTS