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.

My First Salvus Inversion

This micro-tutorial demonstrates how to set up a first full waveform inversion with Salvus.
The simulations are very light and do not require any serious resources.
Some familiary with Python as well as knowledge of how to set up forward simulations with Salvus is assumed.
For an in-depth tutorial of running full waveform inversions with Salvus, please have a look at one of our other tutorials.
Copy
import os
import numpy as np
import xarray as xr

import salvus.namespace as sn

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

Sources and Receivers

We'll generate a ring of transducers, five of which act as sources. This mimics an acquisition setup common in medical devices.
# Generate a ring with 5 scalar sources.
sources = sn.simple_config.source.cartesian.collections.ScalarPoint2DRing(
    x=0, y=0, radius=0.09, count=5, f=1.0
)

events = []
# Use the same receiver ring for each but mask out the receivers close to each
# source.
for i, s in enumerate(sources._sources):
    all_receivers = (
        sn.simple_config.receiver.cartesian.collections.RingPoint2D(
            x=0, y=0, radius=0.09, count=100, fields=["phi"]
        )
    )
    receivers = [
        r
        for r in all_receivers._receivers
        if np.sqrt(
            (s.location[0] - r.location[0]) ** 2
            + (s.location[1] - r.location[1]) ** 2
        )
        > 0.03
    ]
    events.append(
        sn.Event(event_name=f"event_{i}", sources=s, receivers=receivers)
    )
For FWI we'll need something to invert for. For most applications this will be recorded data from some real world setting. Here we'll define a synthetic inversion target and generate our own data.
def get_spherical_inclusion() -> xr.Dataset:
    nx, ny = 200, 200
    x = np.linspace(-0.1, +0.1, nx)
    y = np.linspace(-0.1, +0.1, ny)
    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(
        data_vars={
            "vp": (["x", "y"], vp),
            "rho": (["x", "y"], rho),
        },
        coords={"x": x, "y": y},
    )

    return ds


true_model = get_spherical_inclusion()
ax = true_model.vp.T.plot(figsize=(10, 6))
ax.axes.set_aspect("equal")
# Use the volume model to set up the project.
volume_model = sn.model.volume.cartesian.GenericModel(
    name="true_model", data=true_model
)
p = sn.Project.from_volume_model(
    path="project_full_waveform_inversion",
    volume_model=volume_model,
    load_if_exists=True,
)

# Add the events.
for event in events:
    p.add_to_project(event)

p.visualizations.notebook.domain()
# Generate the data for the target model.
wsc = sn.WaveformSimulationConfiguration(end_time_in_seconds=0.00015)
wsc.physics.wave_equation.time_step_in_seconds = 4.0e-7

ec = sn.EventConfiguration(
    waveform_simulation_configuration=wsc,
    wavelet=sn.simple_config.stf.Ricker(center_frequency=50000.0),
)
p += sn.SimulationConfiguration(
    name="true_model_100kHz",
    elements_per_wavelength=2,
    tensor_order=4,
    max_frequency_in_hertz=100000.0,
    model_configuration=sn.ModelConfiguration(
        background_model=None, volume_models="true_model"
    ),
    event_configuration=ec,
)
p.simulations.launch(
    simulation_configuration="true_model_100kHz",
    events=p.events.list(),
    site_name=SALVUS_FLOW_SITE_NAME,
    ranks_per_job=1,
)
p.simulations.query(block=True)
[2024-05-28 12:38:22,406] INFO: Creating mesh. Hang on.
[2024-05-28 12:38:22,574] INFO: Submitting job array with 5 jobs ...

True
Now we do exactly the same thing, but we'll use a homogeneous initial model.
bm = sn.model.background.homogeneous.IsotropicAcoustic(vp=1500.0, rho=980.0)
mc = sn.ModelConfiguration(background_model=bm)
wsc = sn.WaveformSimulationConfiguration(end_time_in_seconds=0.00015)

p += sn.SimulationConfiguration(
    name="initial_model",
    elements_per_wavelength=2,
    tensor_order=2,
    max_frequency_in_hertz=100000.0,
    model_configuration=mc,
    event_configuration=sn.EventConfiguration(
        waveform_simulation_configuration=wsc,
        wavelet=sn.simple_config.stf.Ricker(center_frequency=50000.0),
    ),
)
p.simulations.launch(
    simulation_configuration="initial_model",
    events=p.events.list(),
    site_name=SALVUS_FLOW_SITE_NAME,
    ranks_per_job=1,
)
p.simulations.query(block=True)
[2024-05-28 12:38:31,818] INFO: Creating mesh. Hang on.
[2024-05-28 12:38:31,902] INFO: Submitting job array with 5 jobs ...

True
Salvus allows for a lot of control and flexibility in how exactly inverse problems are tackled. Things like data processing, window selection, custom weighting schemes, data normalization, bound constraints, wavelength dependent smoothing, and many more are supported but we don't show them in this micro-tutorial. Please refer to some advanced tutorials for more information.
In all cases, an inverse problem requires the definition of a misfit, e.g. the optimization target, as well the actual configuration of the inverse