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
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 proj
if pathlib.Path(PROJECT_DIR).exists():
    print("Opening existing project.")
    p = sn.Project(path=PROJECT_DIR)
else:
    print("Creating new 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)

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

p += e
Creating new project.

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())
[2020-09-28 17:56:14,461] INFO - salvus.project.meshing: Creating mesh. Hang on.
Interpolating model: vs_blob.
/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:702: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  return len(self._mapping)
/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:720: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  yield from self._mapping
When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
<salvus.flow.simple_config.simulation.Waveform object at 0x7fb66b86bad0>
for simulation in ["initial_model", "target_model"]:
    p.simulations.launch(
        ranks_per_job=4,
        site_name=SALVUS_FLOW_SITE_NAME,
        events=p.events.list()[0],
        wall_time_in_seconds_per_job=3600,
        simulation_configuration=simulation,
    )

    while not p.simulations.query():
        time.sleep(1.0)
[2020-09-28 17:56:16,060] INFO - salvus.project.meshing: Creating mesh. Hang on.
[2020-09-28 17:56:16,354] INFO - salvus.project.components.simulation_component: Submitting job ...
Uploading 1 files...

🚀  Submitted [email protected]local
[2020-09-28 17:56:16,880] INFO - salvus.project.components.simulation_component: Status of running jobs (might be slightly out of date).
[2020-09-28 17:56:16,882] INFO - salvus.project.components.simulation_component: job_2009281756476522_977cfbba20:
[2020-09-28 17:56:16,883] INFO - salvus.project.components.simulation_component:     event_SOUTHEAST_OF_LOYALTY_ISLANDS_Mag_5.59_2010-08-03-22-30: pending
[2020-09-28 17:56:17,934] INFO - salvus.project.components.simulation_component: Status of running jobs (might be slightly out of date).
[2020-09-28 17:56:17,935] INFO - salvus.project.components.simulation_component: job_2009281756476522_977cfbba20:
[2020-09-28 17:56:17,936] INFO - salvus.project.components.simulation_component:     event_SOUTHEAST_OF_LOYALTY_ISLANDS_Mag_5.59_2010-08-03-22-30: running `Computing forward wavefield`: 33% done
[2020-09-28 17:56:19,041] INFO - salvus.project.components.simulation_component: Retrieving output for 1 simulation(s) from site.
[2020-09-28 17:56:19,137] INFO - salvus.project.components.simulation_component: Submitting job ...
Uploading 1 files...

🚀  Submitted [email protected]local
[2020-09-28 17:56:19,351] INFO - salvus.project.components.simulation_component: Status of running jobs (might be slightly out of date).
[2020-09-28 17:56:19,352] INFO - salvus.project.components.simulation_component: job_2009281756143736_220874ae0e:
[2020-09-28 17:56:19,353] INFO - salvus.project.components.simulation_component:     event_SOUTHEAST_OF_LOYALTY_ISLANDS_Mag_5.59_2010-08-03-22-30: pending
[2020-09-28 17:56:22,184] INFO - salvus.project.components.simulation_component: Status of running jobs (might be slightly out of date).
[2020-09-28 17:56:22,184] INFO - salvus.project.components.simulation_component: job_2009281756143736_220874ae0e:
[2020-09-28 17:56:22,185] INFO - salvus.project.components.simulation_component:     event_SOUTHEAST_OF_LOYALTY_ISLANDS_Mag_5.59_2010-08-03-22-30: running
[2020-09-28 17:56:23,239] INFO - salvus.project.components.simulation_component: Status of running jobs (might be slightly out of date).
[2020-09-28 17:56:23,240] INFO - salvus.project.components.simulation_component: job_2009281756143736_220874ae0e:
[2020-09-28 17:56:23,241] INFO - salvus.project.components.simulation_component:     event_SOUTHEAST_OF_LOYALTY_ISLANDS_Mag_5.59_2010-08-03-22-30: running `Computing forward wavefield`: 15% done
[2020-09-28 17:56:24,282] INFO - salvus.project.components.simulation_component: Status of running jobs (might be slightly out of date).
[2020-09-28 17:56:24,283] INFO - salvus.project.components.simulation_component: job_2009281756143736_220874ae0e:
[2020-09-28 17:56:24,283] INFO - salvus.project.components.simulation_component:     event_SOUTHEAST_OF_LOYALTY_ISLANDS_Mag_5.59_2010-08-03-22-30: running `Computing forward wavefield`: 80% done
[2020-09-28 17:56:25,389] INFO - salvus.project.components.simulation_component: Status of running jobs (might be slightly out of date).
[2020-09-28 17:56:25,390] INFO - salvus.project.components.simulation_component: job_2009281756143736_220874ae0e:
[2020-09-28 17:56:25,391] INFO - salvus.project.components.simulation_component:     event_SOUTHEAST_OF_LOYALTY_ISLANDS_Mag_5.59_2010-08-03-22-30: finished
[2020-09-28 17:56:26,493] INFO - salvus.project.components.simulation_component: Status of running jobs (might be slightly out of date).
[2020-09-28 17:56:26,494] INFO - salvus.project.components.simulation_component: job_2009281756143736_220874ae0e:
[2020-09-28 17:56:26,495] INFO - salvus.project.components.simulation_component:     event_SOUTHEAST_OF_LOYALTY_ISLANDS_Mag_5.59_2010-08-03-22-30: finished
[2020-09-28 17:56:27,598] INFO - salvus.project.components.simulation_component: Retrieving output for 1 simulation(s) from site.

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,
        },
    )

-

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",
)
When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
p.viz.seismology.windows("initial_model_100_to_300s")
/miniconda/envs/salvus/lib/python3.7/site-packages/ipykernel_launcher.py:1: UserWarning: FixedFormatter should only be used together with FixedLocator
  """Entry point for launching an IPython kernel.
When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.

This examples defines a custom misfit and adjoint source computation but Salvus also already defines a few.

def l2_misfit_and_adjoint_source(
    data_synthetic: np.ndarray,
    data_observed: np.ndarray,
    sampling_rate_in_hertz: float,
):
    # Necessary to not run into crazy small numbers.
    f = 1e10
    assert data_synthetic.shape == data_observed.shape

    # the adjoint source is the negative derivative w.r.t.
    # synthetic data, i.e., we flip the sign here
    diff = data_observed - data_synthetic

    misfit = 0.5 * (diff ** 2).sum() / sampling_rate_in_hertz * f
    adj_src = diff / sampling_rate_in_hertz * f

    return misfit, adj_src


# Add a misfit configuration to the project.
# Note that this is independet of the synthetic data of the forward run
# but it only defines how to compare two sets of data.
p += sn.MisfitConfiguration(
    name="L2-misfit-to-target-model",
    observed_data="PROCESSED_DATA:target_model_with_noise",
    misfit_function=l2_misfit_and_adjoint_source,
    receiver_field="displacement",
    data_selection_configuration="initial_model_100_to_300s",
)

Now we are finally ready to compute misfits and gradients. Note that the misfits can be readily computed without running an additional simulation, whereas the gradients require an upfront adjoint simulation.

# Misfits can be computed before running an adjoint simulation.
misfits = None
while not misfits:
    misfits = p.actions.inversion.compute_misfits(
        simulation_configuration="initial_model",
        misfit_configuration="L2-misfit-to-target-model",
        store_checkpoints=True,
        events=p.events.list()[0],
        ranks_per_job=4,
        site_name=SALVUS_FLOW_SITE_NAME,
    )
    time.sleep(5.0)

print(misfits)
[2020-09-28 17:56:31,753] INFO - salvus.project.components.simulation_component: The following events have been simulated before, but checkpoints are not available for this combination of `site_name` and `ranks_per_job`. They will be run again and the old data will be deleted: ['event_SOUTHEAST_OF_LOYALTY_ISLANDS_Mag_5.59_2010-08-03-22-30']
[2020-09-28 17:56:31,780] INFO - salvus.project.components.simulation_component: Submitting job ...
[2020-09-28 17:56:31,882] INFO - salvus.project.components.action.inversion.inversion_action_component: Launched simulations for 1 events. Please check again to see if they are finished.
[2020-09-28 17:56:37,002] INFO - salvus.project.components.simulation_component: Retrieving output for 1 simulation(s) from site.
{'event_SOUTHEAST_OF_LOYALTY_ISLANDS_Mag_5.59_2010-08-03-22-30': 0.0923226638620134}
# Actually compute the gradients.
while not p.actions.inversion.compute_gradients(
    simulation_configuration="initial_model",
    misfit_configuration="L2-misfit-to-target-model",
    events=p.events.list()[0],
    ranks_per_job=4,
    site_name=SALVUS_FLOW_SITE_NAME,
):
    time.sleep(5.0)
p.viz.nb.gradients(
    simulation_configuration="initial_model",
    misfit_configuration="L2-misfit-to-target-model",
    events=p.events.list(),
)
[2020-09-28 17:56:42,302] INFO - salvus.project.components.simulation_component: Submitting job ...
[2020-09-28 17:56:42,402] INFO - salvus.project.components.action.inversion.inversion_action_component: Launched adjoint simulations for 1 events. Please check again to see if they are finished.
[2020-09-28 17:56:47,537] INFO - salvus.project.components.simulation_component: 1 events have already been submitted. They will not be submitted again.
[2020-09-28 17:56:47,653] INFO - salvus.project.components.action.inversion.inversion_action_component: Some simulations are still running. Please check again to see if they are finished.
[2020-09-28 17:56:53,240] INFO - salvus.project.components.simulation_component: 1 events have already been submitted. They will not be submitted again.
[2020-09-28 17:56:53,354] INFO - salvus.project.components.simulation_component: Retrieving output for 1 simulation(s) from site.
/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:702: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  return len(self._mapping)
/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:720: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  yield from self._mapping
When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.

Salvus guarantees that the gradient is correct with respect to what it is given but what does it actually mean in physical terms? What pieces of data steer and influence the final gradient?

One crucial piece are of course the previously plotted windows and their distribution.

There are a few more ways to analyze the misfits.

# Look at the misfits from the different receivers.
p.viz.nb.misfit_comparison(
    reference_data="initial_model",
    misfit_configuration="L2-misfit-to-target-model",
    event=p.events.list()[0],
)
Misfit comparison for event 'event_SOUTHEAST_OF_LOYALTY_ISLANDS_Mag_5.59_2010-08-03-22-30' and misfit configuration 'L2-misfit-to-target-model'.
initial_model (reference)
IU.CTAO. 4.3717e-03
IU.HNR. 1.1875e-02
IU.MBWA. 2.5338e-02
IU.NWAO. 3.5355e-02
IU.PMG. 5.1214e-03
IU.SNZO. 1.0262e-02
<pandas.io.formats.style.Styler at 0x7fb6650e0a50>
# Can also be shown on a map.
p.viz.seismology.misfit_map(
    reference_data="initial_model",
    misfit_configuration="L2-misfit-to-target-model",
)
When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
# Or have a detailed look at how the individual misfits are computed.
p.viz.nb.misfits(
    simulation_configuration="initial_model",
    misfit_configuration="L2-misfit-to-target-model",
)
When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
p.viz.misfit_histogram(
    simulation_configuration_a="initial_model",
    misfit_configuration="L2-misfit-to-target-model",
    events=p.events.list(),
)
PAGE CONTENTS