Mondaic
PricingAbout UsDocumentation
Mondaic
Version:
No results found.

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.

Lamb's problem

Part 4 - Heterogeneous models

In this part, we consider media with non-constant material parameters and demonstrate how external models can be read into a project.
Copy
%matplotlib inline
# This notebook will use this variable to determine which
# remote site to run on.
import os

SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
PROJECT_DIR = "project"
import numpy as np
import xarray as xr
import salvus.namespace as sn
d = sn.domain.dim2.BoxDomain(x0=0.0, x1=2000.0, y0=0.0, y1=1000.0)

p = sn.Project.from_domain(
    path="project_lambs_problem_volumetric_model",
    domain=d,
    load_if_exists=True,
)

src = sn.simple_config.source.cartesian.MomentTensorPoint2D(
    x=500.0, y=500.0, mxx=1.0, myy=1.0, mxy=0.0
)

recs = [
    sn.simple_config.receiver.cartesian.Point2D(
        y=800.0,
        x=x,
        network_code="REC",
        station_code=f"{_i:05}",
        fields=["displacement"],
    )
    for _i, x in enumerate(np.linspace(1010.0, 1410.0, 5))
]

p.add_to_project(sn.Event(event_name="event_0", sources=src, receivers=recs))
import numpy as np
import xarray as xr

# It helps to define the model larger than your domain.
nx, ny = 600, 400
x = np.linspace(-1000, 3000, nx)
y = np.linspace(-1000, 2000, ny)
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 function.
r = ((xx - 1000.0) / 200) ** 2 + ((yy - 500) / 200) ** 2
vp -= (np.sin(r) / r) * 1000

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

ds.VP.T.plot(figsize=(10, 6))
<matplotlib.collections.QuadMesh at 0x7f0ff03f9650>