This tutorial is not for the latest stable Salvus version (viewing 0.12.16).

%matplotlib inline

import os
import numpy as np
import xarray as xr
import salvus.namespace as sn

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

Simple surface topography

Micro-tutorial

One of the most exciting features of the spectral-element method is the ease with which surface topography or other nonuniform boundaries can be represented in a simulation mesh. In this short tutorial, we will generate a custom topographic model and apply it to a simulation in both 2- and 3-D. The simple trigonometric definition of topography specified below can be expanded to provide much more complexity. Try it out yourself!

2D

Specify domain

d = sn.domain.dim2.BoxDomain(x0=0, x1=20e3, y0=-5e3, y1=0)

Create a project

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

Create a topography model

Generate model

max_x = d.bounding_box[:, 0].max()
x = np.linspace(0, max_x, 10000)
y = np.sin(x * 2 * np.pi / max_x) * 1e3

Wrap in an xarray dataset, specifying the reference elevation

Before application to the mesh, the reference elevation will be subtracted from the topography model in the digital elevation model (DEM).

ds = xr.Dataset(
    data_vars={"dem": (["x"], y, {"reference_elevation": 0.0})},
    coords={"x": x},
)

Turn into a Salvus topography object, plot to validate

tm = sn.topography.cartesian.SurfaceTopography(name="surface", data=ds)
tm.ds.dem.plot()

[<matplotlib.lines.Line2D at 0x7f704663e110>]

Add to your project

p.add_to_project(tm, overwrite=True)

Specify a background model

vp, vs, rho = 4e3, 2.4e3, 2.6e3
bm = sn.model.background.homogeneous.IsotropicElastic(vp=vp, vs=vs, rho=rho)

Create an event for testing

Single source, buried 10m below the surface.

src = sn.simple_config.source.cartesian.SideSetVectorPoint2D(
    point=(5.0e3, 0.0),
    direction="y",
    offset=-10.0,
    side_set_name="y1",
    fx=0.0,
    fy=1.0,
)

11 receivers, located right at the free surface.

recs = [
    sn.simple_config.receiver.cartesian.SideSetPoint2D(
        point=(x, 0.0),
        direction="y",
        offset=0.0,
        side_set_name="y1",
        station_code=f"XX_{i}",
        fields=["velocity", "strain"],
    )
    for i, x in enumerate(np.linspace(12.5e3, 17.5e3, 11))
]

Add to the project

p.add_to_project(sn.EventCollection.from_sources(sources=src, receivers=recs))

Specify a simulation

f_max = 10.0
p.add_to_project(
    sn.SimulationConfiguration(
        tensor_order=2,
        name="sim_with_topo",
        elements_per_wavelength=2.0,
        max_frequency_in_hertz=f_max,
        model_configuration=sn.ModelConfiguration(background_model=bm),
        topography_configuration=sn.TopographyConfiguration("surface"),
        event_configuration=sn.EventConfiguration(
            sn.simple_config.stf.Ricker(center_frequency=f_max / 2),
            sn.WaveformSimulationConfiguration(end_time_in_seconds=10.0),
        ),
        absorbing_boundaries=sn.AbsorbingBoundaryParameters(
            reference_velocity=vp,
            number_of_wavelengths=0.0,
            reference_frequency=f_max / 2,
        ),
    ),
    overwrite=True,
)

Plot the simulation setup

p.viz.nb.simulation_setup("sim_with_topo", events="event_0000")

[2024-03-15 09:05:56,952] INFO: Creating mesh. Hang on.

<salvus.flow.simple_config.simulation.waveform.Waveform object at 0x7f7046694b90>

Run the simulation

p.simulations.launch(
    "sim_with_topo",
    events="event_0000",
    site_name=SALVUS_FLOW_SITE_NAME,
    ranks_per_job=1,
)

p.simulations.query(block=True)

[2024-03-15 09:05:58,070] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2403150905217203_2c053623ba@local

True

Look at some data

p.waveforms.get(
    data_name="SYNTHETIC_DATA:sim_with_topo", events=["event_0000"]
)[0].plot(component="Y", receiver_field="velocity")

3-D

Tutorial coming soon! In the meantime, check out our tutorial on running simulations through a model of Mt St. Helens.