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 1 - First simulation

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

# 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

## Initializing a project and its domain

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 show 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 faciliate 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.
The analytic benchmark we will be comparing the simulated to in part 2 used a 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=-1.0
)
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.
For receivers, we'll set up a line of 5 equidistantly spaced 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 = [
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))
]
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 succesfully 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 now 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.

## Running a simulation

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 influences the wavelengths in the simulation and thus the mesh that discretizes the domain.

#### Defining a model

The analytical solution was computed in an unbounded homogeneous isotropic elastic medium with material parameters specified in SI units as $(\rho = 2200, v_p = 3000, 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) volumetric models for now. So, the next step is to define our background model using the salvus model interface.
mc = sn.ModelConfiguration(
background_model=sn.model.background.homogeneous.IsotropicElastic(
rho=2200.0, vp=3000.0, vs=1847.5
)
)

#### Defining parameters of the source time function

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, as well provide some basic information about the temporal extent of our upcoming simulations. The reference data were computed with using a Ricker wavelet with a center frequency of $15 Hz$ and, looking at the traces plotted above, we can see that the data runs for a bit more than 0.5 seconds. These parameters are now used to define our EventConfiguration object.
ec = sn.EventConfiguration(
wavelet=sn.simple_config.stf.Ricker(center_frequency=15.0),
waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
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 $15Hz$, 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.

#### Defining simulation parameters

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 $30Hz$. 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 $30Hz$. 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,
)
)

#### Visualizing the configuration

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()
)
[2021-07-16 15:30:24,380] INFO: Creating mesh. Hang on. <salvus.flow.simple_config.simulation.Waveform object at 0x7f9be5f0d410>
Feel free to experiement with the dropdown menus and buttons. This visualization can really help debug obvious issues. At this point those of you familiar with older versions of Salvus might be wondering: where did the mesh from? In SalvusProject the complexity of mesh generation is moved into the background, and is handled internally via a reference to the SimulationConfiguration object. While the benefits of this approach are small for small domains and homogeneous models, they will become much later when we consider 3-D models and domains with topography.

#### Running the solver

With everything ready to go, it's now time to run our first simulation! The launch_simulations command below takes a few arguments worth describing:
• site_name: This is an identifier which tells Flow whether you're running on your local machine, some remote cluster, or perhaps the old chess computer in your grandfather's basement. As long as Salvus has been set up correctly on the specified site all data transfers to / from the local or remote machine will happen automatically. Additionally, if a job management system is present on the remote site Flow will monitor the job queue.
• ranks_per_job: This is the number of MPI ranks the job will run on, and can range from 1 to whatever your license will allow.
• events: A list of events for which to run simulations for.
• simulation_configuration: The configuration for which to run simulations for.
p.simulations.launch(
ranks_per_job=2,
site_name=SALVUS_FLOW_SITE_NAME,
events=p.events.list(),
simulation_configuration="my_first_simulation",
)
[2021-07-16 15:30:25,948] INFO: Submitting job ...

1
And that's it! The simulations are off and running. SalvusFlow will take care of abstracting the machine architecture, and SalvusProject will take care of saving all the output data into the correct location, copying it from any remote machines as necessary. We can get the current status of the simulations by calling query_simulations() as below. The parameter block specifies that we wait until all simulations are done.
p.simulations.query(block=True) True
p.viz.nb.waveforms("my_first_simulation", receiver_field="displacement")