Version:

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

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())
[2024-11-06 22:38:39,192] INFO: Creating mesh. Hang on.
Interpolating model: vs_blob.