Version:

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.

Seismological Gradients in SalvusProject

This notebooks explains how to work with gradients with respect to model parameters in a seismological setting inside SalvusProject. We will use a small synthetic example by first creating a model with a velocity perturbation in the center.
Copy
# 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")
PROJECT_DIR = "project"
import time
import pathlib

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

from salvus import namespace as sn
--> Server: 'https://data.mondaic.com/license_server/licensing_server', User: '_MONDAIC_INTERNAL_', Group: '_MONDAIC_INTERNAL_'.
--> Negotiating 1 license instance(s) for 'SalvusMesh' [license version 1.0.0] for 1 seconds ...
--> Success! [Total duration: 0.75 seconds]
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
PROJECT_DIR = "project"

Setup the domain and add an event

# Define a spherical chunk domain.
lat_c, lon_c = -27.5, 135
lat_e, lon_e = 25.5, 42.0

# Uncomment the following line to start from scratch
# !rm -rf project
d = sn.domain.dim3.SphericalChunkDomain(
    lat_center=lat_c,
    lat_extent=lat_e,
    lon_center=lon_c,
    lon_extent=lon_e,
    radius_in_meter=6371e3,
)
p = sn.Project.from_domain(path=PROJECT_DIR, domain=d, load_if_exists=True)

# Add an event by parsing some prepared data.
e = sn.Event(
    sources=sn.simple_config.source.seismology.parse(
        "data/event.txt", dimensions=3
    ),
    receivers=sn.simple_config.receiver.seismology.parse(
        obspy.read_inventory("data/stations.txt"),
        dimensions=3,
        fields=["displacement"],
    ),
)

p += e
The target model has a Gaussian velocity perturbation blob in the middle of. We define it here as an xarray dataset. This can be used to define arbitrary material parameter models that SalvusProject can directly utilize.
lat = np.linspace(lat_c - lat_e / 2.0, lat_c + lat_e / 2.0, 50)
lon = np.linspace(lon_c - lon_e / 2.0, lon_c + lon_e / 2.0, 50)
depth = np.linspace(0, 1000.0, 40)

_, lon_grid, lat_grid = np.meshgrid(depth, lon, lat, indexing="ij")

# Distance from center in degree.
d = ((lon_grid - lon_c) ** 2 + (lat_grid - lat_c) ** 2) ** 0.5

# Apply Gaussian
sigma = 5.0
d = 1.0 / (sigma * np.sqrt(2 * np.pi)) * np.exp(-0.5 * (d / sigma) ** 2)

# Normalize to 10 % deviation
d /= d.max()
d *= 10.0

ds = xr.Dataset(
    data_vars={
        "vs": (["depth", "longitude", "latitude"], d),
        "vp": (["depth", "longitude", "latitude"], d),
    },
    coords={"depth": depth, "latitude": lat, "longitude": lon},
)

ds.vs.attrs["units"] = "%"
ds.vp.attrs["units"] = "%"

# Same for the coordinate axes.
ds.depth.attrs["units"] = "km"
ds.latitude.attrs["units"] = "degrees_north"
ds.longitude.attrs["units"] = "degrees_east"
ds.attrs["geospatial_lat_units"] = "degrees_north"
ds.attrs["geospatial_lon_units"] = "degrees_east"
ds.attrs["geospatial_vertical_units"] = "km"

ds.vs.isel(depth=0).T.plot()
plt.show()
We'll create two simulations: (1) Through the the laterally homogeneous initial model and (2) through the target model with the blob in the middle.
ec = sn.EventConfiguration(
    wavelet=sn.simple_config.stf.GaussianRate(half_duration_in_seconds=100.0),
    waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
        end_time_in_seconds=1800.0
    ),
)
mc = sn.ModelConfiguration(background_model="prem_iso_no_crust")
p += sn.SimulationConfiguration(
    tensor_order=1,
    name="initial_model",
    elements_per_wavelength=1.0,
    min_period_in_seconds=200,
    max_depth_in_meters=500e3,
    model_configuration=sn.ModelConfiguration(
        background_model="prem_iso_no_crust",
    ),
    event_configuration=ec,
)
p += sn.model.volume.seismology.GenericModel(name="vs_blob", data=ds)
p += sn.SimulationConfiguration(
    tensor_order=1,
    name="target_model",
    elements_per_wavelength=1.0,
    min_period_in_seconds=200,
    max_depth_in_meters=500e3,
    model_configuration=sn.ModelConfiguration(
        background_model="prem_iso_no_crust", volume_models=["vs_blob"]
    ),
    event_configuration=ec,
)
p.viz.nb.simulation_setup("target_model", events=p.events.list())
[2021-10-11 10:51:23,912] INFO: Creating mesh. Hang on.
Interpolating model: vs_blob.
<salvus.flow.simple_config.simulation.Waveform object at 0x7f36a59a3b90>
for simulation in ["initial_model", "target_model"]:
    p.simulations.launch(
        ranks_per_job=2,
        site_name=SALVUS_FLOW_SITE_NAME,
        events=p.events.list()[0],
        wall_time_in_seconds_per_job=3600,
        simulation_configuration=simulation,
    )
[2021-10-11 10:51:25,675] INFO: Creating mesh. Hang on.
[2021-10-11 10:51:26,035] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted [email protected]local_cuda
[2021-10-11 10:51:26,346] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted [email protected]local_cuda
p.simulations.query(block=True, ping_interval_in_seconds=1.0)
True
Let's make it a bit more interesting and add some normal distributed noise to the "observed" data through the target mode. The generic processing function approach allows the arbitrary modification of data and it will be applied by SalvusProject on-the-fly.
def process(st, inv, receiver, sources):
    # This particular function adds some random noise to
    # the data. To make it reproducible let's create some
    # deterministic seed - the seed is based on the receiver
    # number and the point-to-point amplitude of the first
    # trace.
    np.random.seed(
        int.from_bytes((receiver.name + str(st[0].data.ptp())).encode(), "big")
        % 2 ** 32
    )
    # Add 5% of normal distributed noise.
    for tr in st:
        tr.data += np.random.normal(
            scale=0.05 * tr.data.ptp(), size=tr.data.shape
        )

    # Filter it again to not have really high-frequencies.
    st.filter("lowpass", freq=1.0 / 80.0, zerophase=True)

    return st


p += sn.processing.seismology.SeismologyProcessingConfiguration(
    name="target_model_with_noise",
    data_source_name="SYNTHETIC_DATA:target_model",
    processing_function=process,
)
# ### Pick Windows
#
# Real world applications usually need some form of data selection. In seismology this oftentimes takes the form of picking windows. Salvus offers some utilities to automatically select these.

if not p.entities.has(
    "data_selection_configuration", "initial_model_100_to_300s"
):
    p.actions.seismology.pick_windows(
        data_selection_configuration="initial_model_100_to_300s",
        observed_data_name="PROCESSED_DATA:target_model_with_noise",
        synthetic_data_name="initial_model",
        events=p.events.list(),
        receiver_field="displacement",
        window_taper_width_in_seconds=50.0,
        window_picking_function="built-in",
        window_picking_function_kwargs={
            "minimum_period_in_seconds": 100.0,
            "maximum_period_in_seconds": 300.0,
        },
    )

Visualize Windows

p.viz.nb.waveforms(
    data=[
        "SYNTHETIC_DATA:initial_model",
        "PROCESSED_DATA:target_model_with_noise",
    ],
    receiver_field="displacement",
    data_selection_configuration="initial_model_100_to_300s",
)