Topography

Jupyter Notebook

Note that everything here is fully automatically tested against - and thus guaranteed to work - only for the latest versions of all the Salvus packages. So please update if anything works differently on your machine.

# Magic
%autosave 0
%matplotlib inline
%config Completer.use_jedi = False

Autosave disabled


Introduction¶

A common use case for spectral-element simulations is for the simulation of wavefields in the presence of surface, subsurface, or seabed topography. In this tutorial we'll * Define a set set of interfaces with topography by using interpolating splines. * Generate a mesh which respects the defined topography. * Set up and run an acoustic simulation through the resulting mesh.

As is the case for many of the other tutorials, this example requires the use of SalvusToolbox, which can be downloaded from here. Let's get started by importing the relevant Python packages.

import matplotlib.pyplot as plt
import numpy as np

import salvus_flow.simple_config as config
import salvus_toolbox.toolbox as st
from salvus_flow import api
from salvus_flow import simple_config as config
from salvus_mesh.optimize_dt import optimize_dt
from salvus_mesh.simple_mesh import rho_from_gardeners, vs_from_poisson


Setting up the model¶

Defining the interpolating splines¶

The first thing we'll do is set up the splines which will define the interfaces in our model. Before we get started though, let's set the x extent of our model, which will be $5$ km.

x_min = 0.0
x_max = 5000.0


To help topography, we'll make use of the get_interpolating_splines function from SalvusToolbox. This function requires a set of interpolation points which can be used to anchor the splines, and these points in turn are defined by a set of locations in $x$ and $y$. In the cell below we define a series of $5$ discontinuities by specifying their interpolation points in each coordiante. We multiply by $1000$ here simply for brevity.

layers_x = [
np.array([0.0, 0.2, 1.0, 2.0, 3.0, 4.0, 4.8, 5.0]) * 1000,
np.array([0.0, 0.2, 1.0, 2.0, 3.0, 4.0, 4.8, 5.0]) * 1000,
np.array([0.0, 0.2, 1.0, 2.0, 3.0, 4.0, 4.8, 5.0]) * 1000,
np.array([0.0, 1.5, 3.5, 5.0]) * 1000,
np.array([0.0, 2.5, 5.0]) * 1000,
]

layers_y = [
np.array([2.0, 2.0, 1.9, 1.7, 2.0, 2.1, 2.0, 2.0]) * 1000,
np.array([1.6, 1.6, 1.5, 1.4, 1.3, 1.4, 1.5, 1.5]) * 1000,
np.array([0.5, 0.5, 0.7, 0.6, 1.1, 0.9, 1.2, 1.2]) * 1000,
np.array([0.2, 0.2, 0.4, 0.4]) * 1000,
np.array([0.0, 0.0, 0.0]) * 1000,
]


We'll also define the p-wave velocities for each layer, and use some helper functions from Salvus Mesh to convert these values to realistic density and s-wave velocities.

# Define p-velocities.
vp = np.array([2000.0, 2500.0, 2800.0, 3200.0])

# Compute vs and rho.
vs = vs_from_poisson(vp)
rho = rho_from_gardeners(vp)


The get_interpolating_splines function also accepts an array of interpolations methods. Internally the function uses interp1d from Scipy, so please check out the documentation here for a complete list of interpolation styles.

interpolation_styles = [
"linear",
"linear",
]


Now we're ready to go ahead and generate our interpolating splines. Given the variables we defined above, we can call the function below to get the horizons of our model.

splines = st.get_interpolating_splines(
layers_x, layers_y, kind=interpolation_styles
)


Below we simply just plot the splines to get an idea of what our model will look like in the end. Indeed the extents and horizons look as we expect, and we can move on to generating the mesh itself. Note that none of the horizons here cross each other, and if this were the case than the subsequent mesh generation would fail with an error. We'll re-visit how to deal with pinched-out layers and crossing horizons in a future tutorial.

# Plot the interfaces.
f = plt.figure(figsize=(10, 5))
x_plot = np.linspace(x_min, x_max)
for top, bot in splines:
plt.plot(x_plot, top(x_plot))
plt.plot(x_plot, bot(x_plot))

plt.xlabel("x (m)")
plt.ylabel("y (m)")
plt.title("Interfaces")

Text(0.5, 1.0, 'Interfaces')


Generating the mesh¶

With our interfaces defined, we can now go ahead and generate the mesh itself. For this case a wrapper is included in the Salvus Toolbox repository. This wrapper takes a series of parameters and attempts to build an optimal mesh with topography, and will also mark each individual region with a flag so it can be identified later. A more detailed description of the arguments one can pass is given below: * x_min, x_max: The minimum and maximum x-values defining the model's horizontal extent. These are the dimensions to be considered before the attachment of absorbing boundaries. If absorbing boundaries are attached the resulting extent of the mesh will grow beyond x_min and x_max, but the solution will not be physical in the extruded regions. * splines: The list of interpolation objects which we have generated above. These could of course be generated independently as well as long as their form is the same. * absorbing_boundaries: A tuple consisting of the side sets which will be made absorbing, along with the number of elements to be extruded in each dimension for the absorbing layer. As previously stated, we recommend an absorbing layer of 3.5+ wavelengths if minimal reflections are desired. * slowest_velocities: A list of the slowest velocities in each layer, ordered the same as the interfaces passed through splines (i.e. top to bottom). These values will be used, along with maximum_frequency and elements_per_wavelength to determine the element sizes in each region. * use_refinements: Create a totally unstrucuted mesh by allowing for vertical refinements. Whether or not this will result in a performance increase is dependent on the model, we recommend try with and without refinements to see what works best for you. * maximum_frequency: The maximum frequency to resolve in the mesh, which will be used in conjunction with elements_per_wavelength.

# Maximum frequency to resolve with elements_per_wavelength.
max_frequency = 20.0

# Generate the mesh
mesh, bnd = st.generate_mesh_from_splines_2d(
x_min=0,
x_max=x_max,
splines=splines,
elements_per_wavelength=2,
maximum_frequency=max_frequency,
use_refinements=True,
slowest_velocities=vs,
absorbing_boundaries=(["x0", "x1", "y0"], 10.0),
)


This function has returned a tuple of two elements. The first entry contains a list of mesh chunks, each of which correspond to a separate region in between two of the defined horizons. The second value returns the minimum size in meters of the extruded boundaries. This value will be used later on when running the actual simulation. We can plot the mesh chunks individually using the build in plot functionality of the mesher. Keep in mind that this may take some time for bigger meshes, so if you are running at higher frequencies you may want to skip of comment out the following cell.

f, axs = plt.subplots(2, 2, figsize=(15, 10), sharex=True)
for _i, (ax, sub_mesh) in enumerate(zip(axs.flatten(), mesh)):

plt.sca(ax)
ax.set_xlabel("x (m)")
ax.set_ylabel("y (m)")
ax.set_title(f"Mesh Chunk {_i + 1}")
sub_mesh.plot(show=False, figure=f, linewidths=0.1)


Our four individual mesh chunks are now properly properly sized and deformed, and are ready to be glued together. Each individual chunk has an elemental variable called region attached to it, and we'll use this variable later on to track which region is which when we're attaching parameters. To sum the mesh regions into one continuous mesh all we need to do is simply call the sum operation on the array of sub-meshes as below...

mesh = np.sum(mesh)


...and now the meshing stage is complete. The cell below will visualize the complete mesh. Note that we have the usual side sets (x0, x1, y0, y1) on the edges of our mesh, but we also have additional side sets which mark the internal discontinuities. These internal side sets can be handy when placing sources and receivers, as we'll explore in a future tutorial.

mesh  # Visualize the mesh.


<salvus_mesh.unstructured_mesh.UnstructuredMesh at 0x7f07b3c08e48>


Attaching parameters¶

Attaching the parameters we defined above is simple because the elements are flagged by region.

nodes = mesh.get_element_nodes()[:, :, 0]
vp_a, vs_a, ro_a = np.ones((3, *nodes.shape))
for _i, (vp_val, vs_val, ro_val) in enumerate(zip(vp, vs, rho)):

# Find which elements are in a given region.
idx = np.where(mesh.elemental_fields["region"] == _i)

# Set parameters in that region to a constant value.
vp_a[idx] = vp_val
vs_a[idx] = vs_val
ro_a[idx] = ro_val

# Attach parameters.
for k, v in zip(["VP", "VS", "RHO"], [vp_a, vs_a, ro_a]):
mesh.attach_field(k, v)

# Attach acoustic / elastic flag.
mesh = st.detect_fluid(mesh)


And of course we can visualize the mesh again to ensure everything is as we expect.

mesh.write_h5("mesh.h5")


Setting up the simulation¶

This section has been convered in other tutorials, so we'll only breifly outline the main steps here. First, create a source object and place it directly on the surface. Note that because the surface is curved this process in non-trivial, and to assist we'll use the SideSetVectorPoint2D class.

# Location of source (this will snap to the closet side-set).
loc = [x_max / 2, 2500.0]

# Ricker wavelet with a center frequency of mf / 2.
stf = config.stf.Ricker(center_frequency=max_frequency / 2)

# Create the source.
source = config.source.cartesian.SideSetVectorPoint2D(
fx=1,
fy=1,
point=loc,
direction="y",
side_set_name="y1",
source_time_function=stf,
)


Now we'll bury an array of 1000 receivers 1m below that same deformed surface. At each receiver we'll save the velocity as well as the measured strain.

receivers = config.receiver.cartesian.SideSetHorizontalPointCollection2D(
x=np.linspace(x_min, x_max, 1000),
offset=-1.0,
side_set_name="y1",
station_code="xx",
fields=["velocity", "strain"],
)


We'll then initialize our Waveform object with the mesh, sources, and receivers we just created. This will also call the source and receiver attachment routines and place those object at their desired side sets and offsets.

w = config.simulation.Waveform(mesh=mesh, sources=source, receivers=receivers)


/miniconda/envs/py37/lib/python3.7/site-packages/salvus_flow/simple_config/simulation/__init__.py:176: UserWarning: If you install the joblib module, the src/rec location process would run in parallel.
cpu_count=cpu_count,


And now we'll set some parameters governing the characteristics of the simulation, validate our chosen inputs, and finally plot the mesh with the sources and receivers to ensure that everything is behaving as expected.

# Set end time of the simulation.
# Start-time will be automatically determined.
w.physics.wave_equation.end_time_in_seconds = 2.0

# Simplified HDF5 point output.
w.output.point_data.format = "hdf5"

# Define coupled Clayton-Enqist / Kosloff
# damping boundaries at the 3 absorbing edges.
absorbing = config.boundary.Absorbing(
width_in_meters=float(bnd),
side_sets=["x0", "y0", "x1"],
taper_amplitude=stf.center_frequency,
)

# Add the boundaries to the parameter file.
w.physics.wave_equation.boundaries = [absorbing]

# Add movie output if you feel inclined
# (beware file size restrictions).
# w.output.volume_data.format = "hdf5"
# w.output.volume_data.fields = ["displacement"]
# w.output.volume_data.filename = "output.h5"
# w.output.volume_data.sampling_interval_in_time_steps = 100

# Ensure salvus_compute will accept our paramters.
w.validate()

# Visualize.
# w


Running the simulation¶

To run the simulation we just need to call Salvus Flow's api.run function and our simulation will run on any site where Salvus has been set up. If this site is remote all data will be copied to and fro as required.

ranks = 2
api.run(
ranks=ranks,
input_file=w,
get_all=True,
overwrite=True,
site_name="local",
output_folder="output",
)

Job job_1907171757933884_3adc09dfd0 running on local with 2 rank(s).
Site information:
* Salvus version: 0.9.8-70-gff9c67d3
* Floating point size: 32


* Downloaded 141.7 MB of results to output.
* Total run time: 56.78 seconds.
* Pure simulation time: 55.76 seconds.

/miniconda/envs/py37/lib/python3.7/site-packages/ipykernel_launcher.py:8: SalvusFlowWarning: This should not happen. Please help the developers to debug this case by sending them this warning message. Don't worry - for you everything should still be fine. Failed to delete folder /root/SalvusFlow/run/job_1907171757933884_3adc09dfd0 due to: [Errno 39] Directory not empty: '/root/SalvusFlow/run/job_1907171757933884_3adc09dfd0'; results of listdir(): SUCCESS; Will now run rm -rf "/root/SalvusFlow/run/job_1907171757933884_3adc09dfd0" to still delete the contents of this folder.

<salvus_flow.sites.salvus_job.SalvusJob at 0x7f07a7e81da0>


Analysing the results¶

Once the simulation is done (2D elastic, high accuracy up to 20 Hz -- should take ~ 1 minute on 2 cores), we can then plot shotgathers of all the fields we've output (vector velocity, as well as the elastic strain).

# Setup up figure.
f = plt.figure(figsize=(15, 10))

# Make plot a bit prettier
plt.setp(ax1.get_yticklabels(), visible=False)
plt.setp(ax3.get_yticklabels(), visible=False)
plt.setp(ax4.get_yticklabels(), visible=False)

# Generate a shotgather from and HDF5 receiver file.

# Use the first call to get the plotting dimensions
vx, dt, extent = st.get_shotgather(rec_file, field="velocity", cmp=0)

# Get the rest of the components we've output
vy, _, _ = st.get_shotgather(rec_file, field="velocity", cmp=1)
sxx, _, _ = st.get_shotgather(rec_file, field="strain", cmp=0)
syy, _, _ = st.get_shotgather(rec_file, field="strain", cmp=1)
sxy, _, _ = st.get_shotgather(rec_file, field="strain", cmp=2)

# Normalize and plot the shotgather.
cv_min, cv_max = 0.05 * vx.min(), 0.05 * vx.max()
sxx_min, sxx_max = 0.01 * sxx.min(), 0.01 * sxx.max()
syy_min, syy_max = 0.01 * syy.min(), 0.01 * syy.max()
sxy_min, sxy_max = 0.01 * sxy.min(), 0.01 * sxy.max()

# Plot the different fields
ax0.imshow(vx, vmin=cv_min, vmax=cv_max, extent=extent, aspect="auto")
ax1.imshow(vy, vmin=cv_min, vmax=cv_max, extent=extent, aspect="auto")
ax2.imshow(sxx, vmin=sxx_min, vmax=sxx_max, extent=extent, aspect="auto")
ax3.imshow(syy, vmin=sxx_min, vmax=sxx_max, extent=extent, aspect="auto")
ax4.imshow(sxy, vmin=sxx_min, vmax=sxx_max, extent=extent, aspect="auto")

# Label plots.
ax0.set_title("$v_x$ (m/s)")
ax1.set_title("$v_y$ (m/s)")
ax2.set_title(r"$\varepsilon_{xx}$")
ax3.set_title(r"$\varepsilon_{yy}$")
ax4.set_title(r"$\varepsilon_{xy}$")

# Label axes.
ax0.set_ylabel("Time (s)")
ax2.set_ylabel("Time (s)")
for a in [ax0, ax1, ax2, ax3, ax4]:
a.set_xlabel("x (m)")