This documentation is not for the latest stable Salvus version.

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.
%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")

Waveform Inversion

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

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

import salvus.namespace as sn

Step 1: Generate the domain and ground truth

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 acquisition. The model for density (RHO) and speed of sound (VP) is created from a structured grid.

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 = 980.0 * np.ones_like(xx)
    mask = np.sqrt(xx ** 2 + yy ** 2) < 0.05
    vp[mask] = 1480.0
    rho[mask] = 1000.0

    mask = np.sqrt(xx ** 2 + (yy - 0.025) ** 2) < 0.015
    vp[mask] = 1550.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(
            "vp": (["x", "y"], vp),
            "rho": (["x", "y"], rho),
        coords={"x": x, "y": y},

    return ds

Let's take a look at the model.

true_model = get_spherical_inclusion()

# Plot the xarray dataset.
plt.figure(figsize=(16, 6))
plt.suptitle("Model with spherical inclusions")

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.

# Uncomment the following line to delete a
# potentially existing project for a fresh start
# !rm -rf project
if pathlib.Path("project").exists():
    print("Opening existing project.")
    p = sn.Project(path="project")
    print("Creating new project.")
    vm = sn.model.volume.cartesian.GenericModel(
        name="true_model", data=true_model
    p = sn.Project.from_volume_model(path="project", volume_model=vm)
Creating new project.

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

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 = [
        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 += 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. Here, we also fix the time step to 400 ns. While this is not strictly necessary and could be automatically detected, of course, it will simplify the data comparison across different simulations.

wsc = sn.WaveformSimulationConfiguration(end_time_in_seconds=0.00015)
wsc.physics.wave_equation.time_step_in_seconds = 4.0e-7

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.

ec = sn.EventConfiguration(
p += sn.SimulationConfiguration(
    # Settings that influence the mesh.
        background_model=None, volume_models="true_model"
    # Potentially event dependent settings.

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

[2020-12-03 18:56:15,493] INFO: Creating mesh. Hang on.
[2020-12-03 18:56:15,953] INFO: Submitting job array with 5 jobs ...
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fa0f0752950>

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.

true_data = p.waveforms.get(

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

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

Alright, enough of these simulations. Now, it is actually time to focus on the inversion.

How do we start? Well, we need a model of course.

What do we already know about the medium? Let's pretend we haven't seen the figures above and have no idea about the phantom. The best prior knowledge we have is that the domain is filled with water and we assume that the speed of sound is 1500 m/s and the density is 980 kg/m^3.

We create a homogeneous background model and set up a simulation configuration that contains the same events as the true data.

bm = sn.model.background.homogeneous.IsotropicAcoustic(vp=1500.0, rho=980.0)
mc = sn.ModelConfiguration(background_model=bm)

p += sn.SimulationConfiguration(
    # Settings that influence the mesh.
    # Potentially event dependent settings.

With the model in hand, we could start iterating right away if we were brave enough. However, in all relevant cases we would want to do some quality checks first before burning many CPU hours. So let's start by just looking at the synthetic waveforms that the initial model produces.

[2020-12-03 18:56:23,097] INFO: Creating mesh. Hang on.
[2020-12-03 18:56:23,175] INFO: Submitting job array with 5 jobs ...