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")
# 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.
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f1a7e977050>
WaveformSimulationobject. This object encapsulates the parameters that we sent to Salvus Compute, and allows us to validate most of the parameters we pass before we begin a simulation. This helps to avoid any typos creeping into our final paramter files. We'll first need to finalize the definition of our source and define our receiver array.
MomentTensor) the simulation would give an error explaining that the given source type could not be attached to any element. We'll also reuse the
stfobject that we generated above for the purposes of mesh design. For more information on different source types, please see our documentation.
source = config.source.cartesian.ScalarPoint2D( source_time_function=wavelet, x=8500.0, y=3490.0, f=1 )
SideSetHorizontalPointCollection2Dreceiver type. This allows us to choose a side-set (here,
"y1"), and to place an array of receivers a certain distance from this side-set. In this case we choose to place the receivers 10 m below the ocean bottom. We also need to specify an array of dynamic fields to record. Here we choose
"gradient-of-phi", which can be related to pressure and particle velocity as outlined here. More information on receiver types can be found in our documentation.
receivers = config.receiver.cartesian.SideSetHorizontalPointCollection2D( x=np.linspace(0.0, 17000.0, 1000), offset=-10.0, station_code="xx", side_set_name="y1", fields=["phi", "gradient-of-phi"], )
config.simulation.Waveformclass. This class allows us to choose and validate many of the relevant simulation parameters. We'll first initialize the object by passing our mesh, along with the sources and receivers we defined above. You'll notice that when executing the following cell, a small progress bar appears and a set of "items" are "located". Here what is happening is that the receivers defined above are being placed exactly 10m away from the side-set
"y1". Given that our domain is rectilinear and has no topography, this process is trivial, and we could of course have simply placed the receivers ourselves as we did for the source. However, once we add topography or bathmetry to our domain (to be outlined in a future tutorial), the automatic point placement can become essential for accurate source and receiver localization. For more information on attaching sources and receivers to deformed meshes, please see our documentation.
s = config.simulation.Waveform(mesh=mesh, sources=source, receivers=receivers)
end_time_in_secondsfields, which are part of the
physics.wave_equationparameter group. If they are not set, Salvus will choose sensible defaults for all three of these parameters. The start and end times will default to the values needed to smoothly simulate the entire source wavelet, and the time step will be computed inside Salvus as per the CFL criterion. Here we set the end time to 8 seconds, and keep the rest of the parameters at their default values.
# Leave start time and time step as their auto-detected defaults. # s.physics.wave_equation.start_time_in_seconds = ? # s.physics.wave_equation.time_step_in_seconds = ? # Set end time. s.physics.wave_equation.end_time_in_seconds = 8.0
HomogeneousDirichletcondition for the top
"y1"boundary. This corresponds to a pressure-free (or free-surface) condition in acoustic media. We'll also now use the absorbing boundary parameters we set above to mock an infinite domain extending from the other side-sets. As mentioned above, we've empirically found good results when the boundaries extend 7 or more elements from the edge of the physical domain, and the taper amplitude is equal to the center frequency of the chosen wavelet.
# Define Dirichlet boundary at top surface (pressure-free). dirichlet = config.boundary.HomogeneousDirichlet(side_sets=["y1"]) # Define coupled Clayton-Enqist / Kosloff # damping boundaries at the the remaining side-sets. absorbing = config.boundary.Absorbing( width_in_meters=685.333, side_sets=["x0", "y0", "x1"], taper_amplitude=wavelet.center_frequency, ) # Add the boundaries to the parameter file. s.physics.wave_equation.boundaries = [absorbing, dirichlet]
s.output.point_data.format = "hdf5"
# s.output.volume_data.format = "hdf5" # s.output.volume_data.filename = "output.h5" # s.output.volume_data.fields = ["phi", "gradient-of-phi"] # s.output.volume_data.sampling_interval_in_time_steps = 100
.validate()on the simulation object. This will ensure that there are no typos in your parameter definitions, and that Salvus will accept the parameters and begin to run. It cannot catch all the possible errors which could occur (for example, if a source or receiver is placed outside of the domain, or if the simulation explodes due to a too-large time step), but it does help catch many common things before you set off a big job. For example, if you one of the fields in the
volume_data.fieldssetting, the validator below would catch this.
api.run( input_file=s, ranks=2, site_name=SALVUS_FLOW_SITE_NAME, get_all=True, output_folder="output", overwrite=False, )
Job `job_2105271711686447_3e37e6f16f` running on `local` with 2 rank(s). Site information: * Salvus version: 0.11.33 * Floating point size: 32
* Downloaded 66.6 MB of results to `output`. * Total run time: 52.82 seconds. * Pure simulation time: 51.97 seconds.
<salvus.flow.sites.salvus_job.SalvusJob at 0x7f1a7cb56a10>
outputin your working directory. All of the files generated during the simulation copied there. If you chose to output volumetric data, you can simply open the
.xdmffile in this output directory to visualize the wavefield in Paraview. Otherwise, you can use the included
salvus.toolboxmodule to generate a shotgather directly from the output file as shown below.
# Generate a shotgather from and HDF5 receiver file. data, dt, extent = toolbox.get_shotgather("output/receivers.h5", field="phi") # Normalize and plot the shotgather. clip_min = 0.01 * data.min() clip_max = 0.01 * data.max() f = plt.figure(figsize=(15, 5)) plt.imshow(data, vmin=clip_min, vmax=clip_max, extent=extent, aspect="auto") plt.xlabel("Offset (m)") plt.ylabel("Time (s)") plt.show()