Mondaic
PricingAbout UsDocumentation
Mondaic
Version:
No results found.

This documentation is not for the latest stable Salvus version.

Lamb's problem

Part 2 - Verification and Benchmark

To validate the solutions generated with Salvus, we will compare our results with semi-analytical ones computed using EX2DDIR.
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 pathlib
import numpy as np
import salvus.namespace as sn
The main purpose of SalvusProject is to manage all data associated with simulations and/or inversions. This ensures reproducibility and avoids computing the same task more than once. The cell below reloads the project state after the first part of the tutorial.
print("Opening existing project.")
p = sn.Project(path=PROJECT_DIR)
Opening existing project.
To verify that we indeed still have all the information available, we can do a few sanity checks like plotting the domain and querying a list of events and simulations.
p.viz.nb.domain()
p.events.list()
['event_0']
p.simulations.list()
['my_first_simulation']
Indeed, all entities from part 1 are still available in the project :)
With our domain initialized and our project reloaded, we can now go right ahead and start preparing some scientific data. We want to compare the simulated data to a semi-analytic solution to Lamb's problem, as described in the introduction of part 1. These data are stored in an HDF5 file named reference_data.h5 in the data directory.
Some data formats, such as ASDF or SEGY describe their data with associated headers. We'll see how to add these types of data in a later tutorial, but in this case we are just reading in raw waveform traces with little to no meta information. Because of this we'll need to assist Salvus a little and tell the project to what events this raw data refers to.
p.waveforms.add_external(
    data_name="reference",
    event="event_0",
    data_filename="data/reference_data.h5",
)
We've already used the function p.viz.nb.waveforms() to visualize seismograms in part 1. To compare the simulated seismograms to the reference data, we just need to pass a list instead of a single simulation name.
p.viz.nb.waveforms(
    ["EXTERNAL_DATA:reference", "my_first_simulation"],
    receiver_field="displacement",
)
Hmm. While the basic shape and behaviour of the waveforms looks good, there are certainly some discrepancies between what we've computed and the reference solution. Now our goal is to see if we can get the signals to match perfectly within the time which they overlap.

Challenge 1: Why doesn't it match?

The most obvious issue with the signals above is that they differ greatly at later times. There's quite an obvious issue here: the reference solution was computed assuming an infinite domain, but the domain we defined is finite in extent. In fact this is the main issue: we need to add absorbing boundaries.
To preserve the stability of the wavefield solution in the presence of complex or anisotropic media, Salvus employs a two-stage approach to absorbing boundaries. First, we apply absorbing boundary conditions at the edge of the mesh as outlined here. These conditions provide good absorbing characteristics for wave impacting the boundary at close to normal incidence, and are sufficient for most cases. If a more substantial absorbing profile is desired, one can also pad the simulated domain with a damping layer. This approach follows that given in this paper. Adding damping layers are advantageous in that they can almost completely cancel any boundary reflections, but do require one to enlarge the computational domain and therefore increase the cost of the resultant simulations. We have found that damping layers provide a good quality / performance tradeoff when 3.5 or more wavelengths are present in the absorbing layer.
In previous versions of Salvus absorbing boundary attachment was unfortunately a manual and tedious process. Fortunately, we now provide an interface to automatically extend the domain in a more user-friendly manner. To activate this feature, we first need to set a few parameters to tell the simulation that we do indeed want a layer of extended absorbing boundaries.
abp = sn.AbsorbingBoundaryParameters(
    reference_velocity=3000.0,
    number_of_wavelengths=3.5,
    reference_frequency=15.0,
)
And now we just simply add a new simulation configuration with the boundaries to our project.
sc = p.entities.get(
    entity_type="simulation_configuration", entity_name="my_first_simulation"
)
p += sn.SimulationConfiguration(
    name="simulation_2",
    max_frequency_in_hertz=sc.max_frequency_in_hertz,
    elements_per_wavelength=sc.elements_per_wavelength,
    model_configuration=sc.model_configuration,
    event_configuration=sc.event_configuration,
    absorbing_boundaries=abp,
)
Note that you can either use += or add_to_project to add new entities to the project.
Visualizing this new configuration, we can see that the mesh has been padded with absorbing boundaries on all sides but the free surface.
p.visualizations.nb.simulation_setup(
    simulation_configuration="simulation_2",
    events=p.events.list(),
)
[2024-03-15 09:01:29,812] INFO: Creating mesh. Hang on.