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.

Manual Gradient Computation

This notebook teaches how to compute gradients with respect to a mesh's material parameters using only SalvusFlow and not SalvusProject. SalvusProject internally handles many of the complexities but the approach presented here grants a bit more flexibility.
# Variable used in the notebook to determine which site
# is used to run the simulations.
import os
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
%matplotlib inline

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

import salvus.namespace as sn

Step 1: Generate target and initial model

This example is purely synthetic so we'll generate a target model with a Gaussian blob perturbation in vsv_s in the center of the domain. This will be used to generate the "observed data" in this case. The initial model (the synthetics for cause of this tutorial) is a homogeneous model without the blob.
# Use xarray to define the blob in vs.
def get_target_model():
    x = np.linspace(-500, 500, 200)
    y = np.linspace(-500, 500, 200)
    xx, yy = np.meshgrid(x, y, indexing="ij")

    # Simple gaussian blob around the center of the domain.
    d = (xx**2 + yy**2) ** 0.5
    vs = np.exp(-0.5 * (d / 100.0) ** 2) * 250.0 + 1000.0

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

    return ds

target_model = get_target_model()
<matplotlib.collections.QuadMesh at 0x7f0ea9593810>
Create a mesh - this is directly used as the initial model.
# Build mesh for the initial model.
initial_model_mesh = sn.simple_mesh.CartesianHomogeneousIsotropicElastic2D(

# Shift.
initial_model_mesh.points[:, 0] +=
initial_model_mesh.points[:, 1] +=
Now interpolate the blob onto it thus creating the target mesh. Have a look at the perturbation in vsv_s.
target_model_mesh = initial_model_mesh.copy()

target_model_mesh = sn.toolbox.interpolate_cartesian(
    data=target_model, mesh=target_model_mesh, extrapolate=True
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f0eb5a30d10>
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.
# Explosive source in one corner of the mesh.
src = sn.simple_config.source.cartesian.MomentTensorPoint2D(

# A centered ring of receivers.
recs = sn.simple_config.receiver.cartesian.collections.RingPoint2D(
    x=0, y=0, radius=200.0, count=50, fields=["displacement"]

# Create two simulation objects.
# (1) The target will serve as observed data.
w_target = sn.simple_config.simulation.Waveform(
    mesh=target_model_mesh, sources=src, receivers=recs
w_target.physics.wave_equation.end_time_in_seconds = 1.0

# (2) The initial model will be used to compute synthetics.
w_initial = sn.simple_config.simulation.Waveform(
    mesh=initial_model_mesh, sources=src, receivers=recs
w_initial.physics.wave_equation.end_time_in_seconds = 1.0

# Must store the checkpoints for the subsequent adjoint run.
w_initial.output.volume_data.format = "hdf5"
w_initial.output.volume_data.filename = "output.h5"
w_initial.output.volume_data.fields = ["adjoint-checkpoint"]
w_initial.output.volume_data.sampling_interval_in_time_steps = (
For the target/observed model and the initial/current model.
# Create the "observed data".
j_target =

# Forward run for the synthetics.
j_initial =
    # Don't delete the remote files as they are needed
    # for the adjoint run!
SalvusJob `job_2407051019737769_c05d1130c1` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 2024.1.1
  * Floating point size: 32
-> Current Task: Time loop complete* Downloaded 341.0 KB of results to `output_target`.
* Total run time: 0.81 seconds.
* Pure simulation time: 0.42 seconds.
SalvusJob `job_2407051019399236_5f077f3315` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 2024.1.1
  * Floating point size: 32
-> Current Task: Time loop complete* Downloaded 341.4 KB of results to `output_initial`.
* Total run time: 0.63 seconds.
* Pure simulation time: 0.45 seconds.
# Load the data from both runs, and visualize.
# Note that this is a vectorial receiver and thus the component
# has to be specified.
target_event = j_target.get_as_event()
initial_event = j_initial.get_as_event()

target_event.plot(receiver_field="displacement", component="X")
initial_event.plot(receiver_field="displacement", component="X")