This documentation is not for the latest stable Salvus version.
Pythonpackages into our notebook.
%matplotlib inline %config Completer.use_jedi = False from pathlib import Path import os import requests import shutil from typing import List import scipy import obspy import numpy as np import xarray as xr import matplotlib.pyplot as plt from salvus.flow import api from salvus.flow import simple_config as config from salvus.toolbox import toolbox SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
--> Server: 'https://data.mondaic.com/license_server/licensing_server', User: '_MONDAIC_INTERNAL_', Group: '_MONDAIC_INTERNAL_'. --> Negotiating 1 license instance(s) for 'SalvusMesh' [license version 1.0.0] for 1 seconds ... --> Success! [Total duration: 0.70 seconds]
# Download Marmousi model ~150 MB. target_path = "elastic-marmousi-model.tar.gz" if not Path(target_path).exists(): url = ( "https://s3.amazonaws.com/open.source.geoscience/open_data" "/elastic-marmousi/elastic-marmousi-model.tar.gz" ) response = requests.get(url, stream=True) if response.status_code == 200: with open(target_path, "wb") as f: f.write(response.raw.read()) shutil.unpack_archive(target_path)
xarrayas an intermediary to encapsulate both the model parameters and geometry. More info on
xarray, including extensive documentation, can be found here. Reading our Marmousi model into an
xarray.Datasetis trivial, and the process can be inspected by opening the acompanying
xarraydataset with the following information
# Read marmousi model. model_directory = Path("./elastic-marmousi-model/model") marmousi_model = toolbox.read_elastic_marmousi(model_directory, ["VP", "RHO"])
xarray.Datasetobject is the ability to quickly plot the contained data. Let's do this below as a way to quickly QC that the model reading went correctly.
marmousi_model["VP"].T.plot(aspect="auto", figsize=(15, 5)) plt.show() marmousi_model["RHO"].T.plot(aspect="auto", figsize=(15, 5)) plt.show()
wavelet = config.stf.Ricker(center_frequency=5.0) wavelet.plot()
num_absorbing_layers = 10 absorbing_side_sets = ["x0", "x1", "y0"]
mesh_frequencyparameter to twice the center frequency of the Ricker wavelet we will use, as described above.
mesh_frequency = 2 * wavelet.center_frequency mesh = toolbox.mesh_from_xarray( model_order=4, data=marmousi_model, slowest_velocity="VP", maximum_frequency=mesh_frequency, elements_per_wavelength=1.5, absorbing_boundaries=(absorbing_side_sets, num_absorbing_layers), )
mesh, it's now enough to just type
meshinto the cell below. When run in a notebook, this command will plot an interactive 2- or 3-D plot of the simulation mesh, with the available parameters interpolated to the GLL points defining the model interpolation order. This is why we see variability within each element which is plotted below. You can expand the
Parameterdialog box in order to visualize any parameter present in the mesh. Obviously, we have the spatially variable
RHO, which define 2-D acoustic parameters of our medium. In addition, we have the elemental-parameter
fluid, along with
fluidflag is read by Salvus Compute and is used to determine which elements in the simulation are acoustic and which are elastic. In this case we see that
fluid = 1everywhere in the mesh, which is telling us that the entire mesh should simulation the fluid (acoustic) wave equation. In this tutorial and this tutorial we explore how to run coupled simulations, and see how one can change the fluid flag at will, and even vary it randomly between individual elements. The
absorbing_boundariesflag simply shows us where the domain has been un-physically extended for use by the damping layers. This is handy to know, as we don't want to put sources or receivers in this area. As we specified
["x0", "x1", "y0"]as our absorbing side-sets, you should see the damping layers extending out of the left, bottom, and right sides of the mesh. Since we will be simulating a free-surface condition at the top of the model, no extension has been made here. As a reminder from earlier tutorials: you can visualize the side sets in the mesh by clicking the "Show Side Sets" button in the widget below.