Mondaic - Full Waveform Solutions

An accurate solution to the wave equation is a requirement for a wide variety of seismological research. In this tutorial, we will validate the accuracy of Salvus by comparing numerically calculated seismograms to semi-analytical solutions of Lamb's Problem in 2-D. In addition to giving us confidence in the synthetic data we will use in future tutorials, it also gives us a chance to gently learn some of the key features of the SalvusProject API.

Lamb's problem is concerned with the behavior of the elastic wave equation in the presence of a half-space bounded by a free-surface condition. In our solution we expect both direct arrivals and those reflected from the free-surface, along with a contribution from the 2-D Rayleigh wave. We'll consider a half-space bounded at

y=2000

, and excite waves using a Ricker source with a center frequency of 15 Hz. This setup keeps compute times very low, while also allowing for a fair amount of wavelengths to propagate within our domain. To get started, let's first import all the Python tools we'll need.

%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"

# Uncomment the following line to delete a
# potentially existing project for a fresh start
# !rm -rf project

import pathlib
import numpy as np
import salvus.namespace as sn

Before we initialize our project, we'll first need to initialize the spatial domain to which our project corresponds. In this case we'll be using a simple 2-D box domain. Don't worry, we'll make things a bit more exciting in future tutorials. The box domain can easily be constructed from a set of two-dimensional extents as shown in the cell below.

d = sn.domain.dim2.BoxDomain(x0=0.0, x1=2000.0, y0=0.0, y1=1000.0)

For a given project the domain we specify is immutable; its extents and characteristics are used to infer other information regarding our meshes, simulations, and data. With the simple domain definition given above, we're now ready to initialize our project. To do this we can use the Project.from_domain() constructor as shown below. This function takes a path and a domain object such as the one we just constructed. Calling this function several times will load the project state from disk instead. The project state is persistent and not tied to this notebook, so you can easily take a break and continue working on the project anytime.

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

Most entities of a project can be visualized in a notebook widget by calling p.viz.nb. followed by the name of the entity. Currently the project does not contain a lot of data yet, but at least we can visualize the domain extent.

p.viz.nb.domain()

Next, we need to set up a source and a few receivers. To facilitate this we'll use several simple_config objects from SalvusFlow above. You can use the tab key to browse through the different categories and to see a list of all available options, or check the documentation here.

In the benchmark in part 2, we will be comparing the simulated seismograms to the semi-analytic response of 2-D vector point source with components

[0, -1e10]

at a position of

[1000, 500]

. Here, forces are specified in Newton-meters, with positions in meters.

src = sn.simple_config.source.cartesian.VectorPoint2D(
    x=1000.0, y=500.0, fx=0.0, fy=-1e10
)

Note that we only specify the spatial characteristics of the source at this time. The temporal parameters, most importantly the frequency content of the source wavelet, is tightly coupled to the discretization of the domain. It is common to run several simulations using different frequency bands, which is why those parameters will be defined later as part of the SimulationConfiguration.

As measurement locations, we'll set up a line of 5 equidistantly spaced receivers between

x=1010.0

m and

x=1410.0

m, at a y position of 800 m. We need to specify a unique station code for each receiver, as well as a list of which fields we'd like the receivers to record (displacement in this case).

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

With our sources and receivers now defined, we can add the combination of them as an Event to our project. Every event requires a unique name.

p.add_to_project(sn.Event(event_name="event_0", sources=src, receivers=recs))

Note here the syntax we used. An Event, along with several other relevant objects, can be added to a project by calling add_to_project or simply using the += operator. Once the object is successfully added to the project it is then "serialized", or saved, within the project directory structure. The power and usefulness of this concept will become apparent in a later tutorial. For now all you need to know is that the event is officially a part of our project!

We can verify this by plotting the domain again.

p.viz.nb.domain()

That's enough setup for now. Let's get going with some simulations.

What is still missing for our first simulation? We have not yet specified the material properties within our domain, nor the time characteristics of the source. Both influence the wavelengths in the simulation and thus the mesh that discretizes the domain.

The analytical solution was computed in an unbounded homogeneous isotropic elastic medium with material parameters specified in SI units as

(\rho = 2200, v_p = 3200, v_s = 1847.5)

. As the analytic solution was computed in a homogeneous medium, we don't need to concern ourselves with (2- or 3-D) heterogeneities of volumetric models for now. So, the next step is to define our homogeneous medium using the Salvus interface for a background_model.

mc = sn.ModelConfiguration(
    background_model=sn.model.background.homogeneous.IsotropicElastic(
        rho=2200.0, vp=3200.0, vs=1847.5
    )
)

Note that up until now we have not specified any information regarding the frequency content of the data we are planning on simulating, and in fact all the parameters we've specified have been frequency-independent. This is deliberate, as it is often the case that information on material parameters are provided independent of frequency. The next step is to add a time-frequency axis to our project, which enters in the form of an EventConfiguration. Here, at a bare minimum, we need to specify what type of source wavelet we would like to model, and to provide some basic information about the temporal extent of our upcoming simulations. The reference data were computed with using a Ricker wavelet centered around time

t=0

and with a center frequency of

14.5

Hz. The time interval spans a bit more than 0.5 seconds. These parameters are now used to define our EventConfiguration object. The start time is optional and only necessary here for comparing with the analytic solution in the next part. Salvus will auto-derive an appropriate start time if not given. Because the Ricker wavelet is centered around zero, we have non-zero contributions before

t=0

and consequently need to start the simulation at a negative starting time.

ec = sn.EventConfiguration(
    wavelet=sn.simple_config.stf.Ricker(center_frequency=14.5),
    waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
        start_time_in_seconds=-0.08,
        end_time_in_seconds=0.6,
    ),
)

To get a better sense of what our wavelet looks like in both the time and frequency domain, we can easily plot its characteristics in the cell below.

ec.wavelet.plot()

We quickly see that, while the center frequency of the wavelet was specified to be

14.5

Hz, there is actually a fair bit of energy that exists at frequencies higher than this. Its important to design our simulations so that they properly resolve all the frequencies we are interested in.

The final step in defining a simulation is pulling together all the above into a single reproducible SimulationConfiguration. A SimulationConfiguration is a unique identifier that brings together the model, the source wavelet parameterization, and a proxy of the resolution of the simulation together. Depending on the application, sufficient accuracy is typically achieved with 7 to 10 grid points per wavelength, sometimes even less. By default, we use a polynomial degree of 4 in the spectral-element mesh, which gives 5 grid points per element. Here, we need to specify elements_per_wavelength and choose 1.0.

The question then remains: given a broadband source wavelet, which frequency do we want to mesh for? The wavelet plot above gives us a clue: the vast majority of the energy in the current wavelet is contained at frequencies below

30

Hz. For our first attempt at matching the analytic solution then, we'll require that our mesh be generated using one element per wavelength at a frequency of

30

Hz. As you are probably becoming familiar with by now, we can add the relevant SimulationConfiguration to our project as below.

p.add_to_project(
    sn.SimulationConfiguration(
        name="my_first_simulation",
        max_frequency_in_hertz=30.0,
        elements_per_wavelength=1.0,
        model_configuration=mc,
        event_configuration=ec,
    )
)

So far, regarding our simulation, we have:

Defined a homogeneous material model
Defined a Ricker wavelet source
Set a resolution criterion

In fact, this is all we need to do! Before we actually run the simulation though, it can be helpful to get a visual overview of what is about to happen. Salvus project provides a small convenience function to visualize a SimulationConfiguration directly in the notebook, as below. This function takes a list of events as well, for the purpose of overplotting sources and receivers on the resultant domain. Let's have a look.

p.viz.nb.simulation_setup(
    simulation_configuration="my_first_simulation", events=p.events.list()
)

[2024-07-05 09:52:17,254] INFO: Creating mesh. Hang on.