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.

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 0x7fb428da5350>

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.unstructured_mesh.UnstructuredMesh at 0x7fb421d36250>

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!
Job `job_2104062136035929_b5d4927ae2` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 0.11.29
  * Floating point size: 32
-> Current Task: Time loop complete* Downloaded 342.0 KB of results to `output_target`.
* Total run time: 0.93 seconds.
* Pure simulation time: 0.64 seconds.
Job `job_2104062136854942_9b91c89450` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 0.11.29
  * Floating point size: 32
-> Current Task: Time loop complete* Downloaded 342.4 KB of results to `output_initial`.
* Total run time: 1.08 seconds.
* Pure simulation time: 0.91 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")