Version:

My First Salvus Simulation

This micro-tutorial demonstrates how to set up a first Salvus simulation and teaches some of the fundamental concepts of Salvus.
The simulations are very light and do not require any serious resources. Some familiarity with Python is assumed.
For an in-depth tutorial of running simulations with Salvus, please have a look at one of our other tutorials.
Copy
# Three standard Python imports
import os

import numpy as np
import xarray as xr

# The Salvus package.
import salvus.namespace as sn

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

SalvusProject

SalvusProject is a file structure on disk that organizes and stores simulations and inversions for a given domain.
Salvus supports domains in two and three dimensions for spherical and cartesian coordinates systems. For this micro-tutorial we choose a two-dimensional cartesian box domain.
domain = sn.domain.dim2.BoxDomain(x0=0.0, x1=2000.0, y0=0.0, y1=1000.0)
domain.plot()
# Once we have a domain we can initialize a project with a path on disk.
p = sn.Project.from_domain(
    path="my_project", domain=domain, load_if_exists=True
)
Sources are necessary to inject energy into waveform simulations and most applications require the wavefield sampled at point receivers. In Salvus, both are specified as small Python objects.
# Source and receivers exist for multiple dimensions and coordinate systems.
source = sn.simple_config.source.cartesian.VectorPoint2D(
    x=600.0, y=700.0, fx=0.0, fy=-1e10
)

receivers = [
    sn.simple_config.receiver.cartesian.Point2D(
        y=800.0,
        x=x,
        station_code=f"REC{i + 1}",
        # Note that one can specify the desired recording field here.
        fields=["displacement"],
    )
    for i, x in enumerate(np.linspace(1000.0, 1400.0, 5))
]

# The source and receivers object are quite similar to Python dictionaries.
source, receivers
({'location': [600.0, 700.0],
  'spatial_type': 'vector',
  'spatial_weights': [0.0, -10000000000.0]},
 [{'fields': ['displacement'],
   'location': [1000.0, 800.0],
   'location_code': '',
   'network_code': 'XX',
   'station_code': 'REC1'},
  {'fields': ['displacement'],
   'location': [1100.0, 800.0],
   'location_code': '',
   'network_code': 'XX',
   'station_code': 'REC2'},
  {'fields': ['displacement'],
   'location': [1200.0, 800.0],
   'location_code': '',
   'network_code': 'XX',
   'station_code': 'REC3'},
  {'fields': ['displacement'],
   'location': [1300.0, 800.0],
   'location_code': '',
   'network_code': 'XX',
   'station_code': 'REC4'},
  {'fields': ['displacement'],
   'location': [1400.0, 800.0],
   'location_code': '',
   'network_code': 'XX',
   'station_code': 'REC5'}])
# Sources and receivers are tied together in an event. An event can contain an
# arbitrary number of sources and receivers and each event will be one simulation.
event = sn.Event(event_name="event_a", sources=source, receivers=receivers)
# Lastly add the event to the already existing project and visualize it.
p.add_to_project(event)
p.visualizations.notebook.domain()
Salvus is very flexible when it comes to specifying material models or numerical meshes and supports models ranging in complexity from homogeneous materials all the way to fully unstructured custom meshes with a variety of wave physics.
To keep it somewhat interesting, we'll specify a heterogeneous 2-D elastic model here. Salvus interfaces with the xarray library to achieve that. This is a widely used scientific Python library and readers/write to a large variety of file formats exists. Also note that nothing in the next cell is Salvus specific code.
# It helps to define the model larger than your domain to avoid
# edge effects during the model interpolation step.
nx = 500
x = np.linspace(-100, 2100, nx)
y = np.linspace(-100, 1100, nx // 2)
xx, yy = np.meshgrid(x, y, indexing="ij")

# Background model.
vp = 3000.0 * np.ones_like(xx)
vs = 1847.5 * np.ones_like(xx)
rho = 2200.0 * np.ones_like(xx)

# Sinc perturbations.
r = ((xx - 1000.0) / 200) ** 2 + ((yy - 500) / 200) ** 2
vp -= (np.sin(r) / r) * 1000

# Assemble into one dataset.
ds = xr.Dataset(
    data_vars={
        "VP": (["x", "y"], vp),
        "VS": (["x", "y"], vs),
        "RHO": (["x", "y"], rho),
    },
    coords={"x": x, "y": y},
)

# Plot it.
ax = ds.VP.T.plot(figsize=(10, 6))
ax.axes.set_aspect("equal")