# Layered Media

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

%matplotlib inline

# Import modules and setup paths.
import os
import h5py
import pyasdf
import matplotlib.pyplot as plt

PARAVIEW_BIN = os.environ["PARAVIEW_BIN"]
SALVUS_BIN = os.environ["SALVUS_BIN"]

all_rec_names = ['rec1', 'rec2', 'rec3', 'rec4', 'rec5']


In this example, we will simulate wave propagation through acoustic, elastic, and coupled media. By the end of the tutorial you should be able to generate a mesh, attach sources and receivers, generate simple movies, and run simulations in layered, coupled fluid / solid media. Let's get started!

### Generating the mesh¶

The first step is to generate a mesh approprate for the problem. These meshes can be generated from commerical meshing packages, such as Trelis, or with the Salvus ecosystem itself, using salvus_meshser. For this tutorial, we'll use Salvus' mesher to generate some simple layered models. Open the file up in your favourite text editor and check out the structure. These .bm mesher files make generating layered models a breeze. You can run the mesher in the box below by typing python -m pymesher.interface. Run it with --help to get information on the options.

In this example, we'll mock a simple oceanic seismic survey. The .bm file specifies a model depth of about 4 km, and let's set the horizontal distance to 10 km. Since we're looking for something that can be run on a laptop, let's keep the mesh relatively small. With those dimensions, a mesh acurate to 10 Hz fulfils this criteria well. Of course, if you've got some beefier hardware, feel free to increase this frequency.

# Edit if necessary!

!python -m pymesher.interface Cartesian2D --basic.model ACOUSTIC_LAYERED.bm \
--basic.period 0.1 --cartesian2D.x 0 10 --overwrite --output_filename acoustic_layered.e

Setting up background model and element sizes...
Creating the skeleton...
Creating the unstructured mesh...
attaching elastic parameters
Computing mesh quality...
Computing resolved period...
Writing mesh to file...
==============================================================================
SUMMARY OF MESH PROPERTIES:

model name                       | ACOUSTIC_LAYERED
dominant period input            |      0.10 s
elements per wavelength          |      2.00
Courant Number                   |      0.60

resolved period (global max)     |      0.10 s
location (z coordinate)        |      1.95 km
resolved period (percentile 95)  |      0.10 s
time step dt                     |    0.0087 s
location (z coordinate)        |      2.55 km
number of elements               |      2412
number of points                 |      2536
cost factor (nelem / dt)         |  2.76e+05

max edge aspect ratio            |      2.68
max equiangular skewness         |      0.59
==============================================================================
GLOBAL VARIABLES:
dt                         |   0.00874
minimum_period             |   0.09976
==============================================================================
ELEMENTAL FIELDS:
RHO_0
RHO_1
RHO_2
RHO_3
VP_0
VP_1
VP_2
VP_3
VS_0
VS_1
VS_2
VS_3
dt
edge_aspect_ratio
element_type
equiangular_skewness
fluid
minimum_period
==============================================================================
NODAL FIELDS:
==============================================================================
SIDE SETS:
x0
x1
y0
y1
==============================================================================

SUCCESSFULLY GENERATED MESH IN 0.359494 SECONDS.
SAVED TO 'acoustic_layered.e' (151.7 KB).


Great. Now let's open it up in Paraview. You should see something similar to the figure below.

!$PARAVIEW_BIN acoustic_layered.e  Now, we'll leave you alone. See if you can generate some seismograms and visualize them here. Remember in this case we're in an acoustic medium! See if you can generate a nice movie, and get some receiver output. Let's run a simulation for 2 seconds. !$SALVUS_BIN \
--mesh-file acoustic_layered.e --model-file acoustic_layered.e \
--polynomial-order 4 --dimension 2 --end-time 2.00 \
--cmd-line-source --source-temporal-type ricker --source-spatial-type scalar \
--source-location-x 6000 --source-location-y 3000 \
--source-scale 1e10 --source-center-frequency 5.0 \
--absorbing-boundaries x0,x1,y0 \
--save-movie --movie-file-name movie.h5 --movie-fields u_ACOUSTIC

====================================================
Salvus version 0.0.1-686-g644f5f44-dirty
Floating point size: 64
Compiled for GLL orders: 4
Running on 1 rank(s).
====================================================

--source-time-delay not set. Will auto-calculate.
Initializing problem.
Start time set to -0.311879.
Time step set to: 0.000748214 s.
Start time set to -0.311879.
Begin time loop.
Time loop progress [100%].
Time loop completed in 12.2517 seconds.
Begin post processing.
Problem complete.

!python petsc_gen_xdmf.py movie.h5
!$PARAVIEW_BIN ./movie.xmf  def plot_hdf5_seismograms(asdf_file, axes, mode): if mode == 'acoustic': par_map = ['0'] elif mode == 'elastic': par_map = ['Z', 'E'] else: raise ValueError("Incorrect mode specified.") flat_axes = axes.flatten() for i, waveform in enumerate(asdf_file.waveforms): for j in range(len(par_map)): axes[i+j*len(par_map)].plot( waveform.displacement.select(component=par_map[j])[0].data) salvus_rec = pyasdf.ASDFDataSet('receivers.h5') fig, axes = plt.subplots(1, 5, sharex='row', sharey='col', figsize=(20,5)) plot_hdf5_seismograms(salvus_rec, axes, 'acoustic')  ### Switching to elastic media¶ Let's up the complexity of our simulation a bit, and add a simulation using elastic media. Use a Poisson's ratio of 0.25 to generate non-zero VS values from your chosen VP values. Write a python function to transform the values in your VP column to VS. # Generating the VS values. def VP_to_VS(VP): """Take a VP value and return VS.""" ratio = 0.25 VS = ((1 - 2 * ratio) / (2 * (1 - ratio))) * VP return VS VS = [VP_to_VS(V) for V in [1.5e3, 1.5e3, 2.7e3, 3.0e4, 3.9e3, 5.0e3, 5.5e3, 6.0e3, 6.2e3]] print("S-wave velocities: ", VS)  S-wave velocities: [500.0, 500.0, 900.0, 10000.0, 1300.0, 1666.6666666666665, 1833.3333333333333, 2000.0, 2066.6666666666665]  Now, use the command line functionality of the mesher to generate an elastic mesh with these new values. As well, let's open it up in Paraview to see what it looks like. # Generating the elastic mesh. # Edit if required. # !mousepad ELASTIC_LAYERED.bm ! python -m pymesher.interface Cartesian2D --basic.model ELASTIC_LAYERED.bm \ --basic.period 0.1 --cartesian2D.x 0 10 \ --overwrite --output_filename elastic_layered.e  Setting up background model and element sizes... Creating the skeleton... Creating the unstructured mesh... attaching elastic parameters Computing mesh quality... Computing resolved period... Writing mesh to file... ============================================================================== SUMMARY OF MESH PROPERTIES: model name | ELASTIC_LAYERED dominant period input | 0.10 s elements per wavelength | 2.00 Courant Number | 0.60 resolved period (global max) | 0.10 s location (z coordinate) | 3.81 km resolved period (percentile 95) | 0.10 s time step dt | 0.0030 s location (z coordinate) | 2.64 km number of elements | 17900 number of points | 18216 cost factor (nelem / dt) | 6.02e+06 max edge aspect ratio | 3.00 max equiangular skewness | 0.63 ============================================================================== GLOBAL VARIABLES: dt | 0.00297 minimum_period | 0.10000 ============================================================================== ELEMENTAL FIELDS: RHO_0 RHO_1 RHO_2 RHO_3 VP_0 VP_1 VP_2 VP_3 VS_0 VS_1 VS_2 VS_3 dt edge_aspect_ratio element_type equiangular_skewness fluid minimum_period ============================================================================== NODAL FIELDS: ============================================================================== SIDE SETS: x0 x1 y0 y1 ============================================================================== SUCCESSFULLY GENERATED MESH IN 0.184591 SECONDS. SAVED TO 'elastic_layered.e' (307.4 KB).  # View in paraview. !$PARAVIEW_BIN ./elastic_layered.e


What do you notice is different about the mesh?

Remember what is happening here. To be accurate at 4th order, Salvus requires a mesh with 2 elements per minimum wavelength. Since the shear-wave velocity is now so much smaller, the elements must be made smaller to compensate. Thus, there are many more elements in this mesh than the equivalent acoustic one. As well since the elements are smaller, the timestep must be reduced to remain stable. All in all an elastic simulation can be hundreds of times more expensive than an equivalent acoustic-only simulation. It is no wonder why most of exploration seismology is done under the acoustic approximation!

If you have decent hardware, a 0.5 second elastic simulation should take just over 3 minutes on one core. Depending on how well your virtualbox supports MPI, you can speed this up by specifying mpirun -n 2. Of course, if we were on a cluster, we could specify mpirun -n N, with N as large as we can support. Alternatively, you can reduce the resolution of the mesh to generate less elements. Don't forget to adjust the center frequency of your source as well!

# Running the elastic simulation
# This is a very short simulation for speed reasons. If you can spare a
# minute or two, make the end-time larger (1.0 seconds or more).
!$SALVUS_BIN \ --mesh-file elastic_layered.e --model-file elastic_layered.e \ --polynomial-order 4 --dimension 2 --end-time 0.01 \ --cmd-line-source \ --source-temporal-type ricker --source-spatial-type moment_tensor \ --source-location-x 6000 --source-location-y 3000 \ --source-scale 1e10,1e10,0 --source-center-frequency 5.0 --number-of-receivers 5 \ --receiver-names rec1,rec2,rec3,rec4,rec5 --receiver-location-x 10,110,210,310,410 \ --receiver-location-y 1800,1800,1800,1800,1800 \ --receiver-file-name receivers.h5 --receiver-fields u_ELASTIC \ --receiver-medium solid,solid,solid,solid,solid --absorbing-boundaries x0,x1,y0 \ --save-movie --movie-file-name movie.h5 --movie-fields u_ELASTIC --save-frame-every 100  ==================================================== Salvus version 0.0.1-686-g644f5f44-dirty Floating point size: 64 Compiled for GLL orders: 4 Running on 1 rank(s). ==================================================== --source-time-delay not set. Will auto-calculate. Initializing problem. Start time set to -0.311879. Time step set to: 0.000255812 s. Start time set to -0.311879. Begin time loop. Time loop progress [100%]. Time loop completed in 56.9428 seconds. Begin post processing. Problem complete.  ### Switching to coupled acoustic / elastic media¶ To really get the physics correct in this example, we'll want to work with a coupled domain. The topmost layer we will treat as an ocean, while the rest of the problem will still be elastic. To generate this mesh, you simply need to specify VS = 0.0 in the top layer of the model file. Make sure you set it to 0.0 at both the top and bottom of the discontinuity, otherwise you will have a gradient in S-wave velocity, and the mesher cannot assume that the element is fluid. If the mesh is generated with fluid and solid elements, Salvus with use PETSc's DMPLEX package to explore the Hasse Diagram for each element, and automatically determine which edges lie along coupling boundaries. No need to specify these manually! # Generating the elastic mesh. # Edit if necessary # !mousepad COUPLED_LAYERED.bm !python -m pymesher.interface Cartesian2D --basic.model COUPLED_LAYERED.bm \ --basic.period 0.1 --cartesian2D.x 0 10 --overwrite --output_filename coupled_layered.e  Setting up background model and element sizes... Creating the skeleton... Creating the unstructured mesh... attaching elastic parameters Computing mesh quality... Computing resolved period... Writing mesh to file... ============================================================================== SUMMARY OF MESH PROPERTIES: model name | COUPLED_LAYERED dominant period input | 0.10 s elements per wavelength | 2.00 Courant Number | 0.60 resolved period (global max) | 0.10 s location (z coordinate) | 2.06 km resolved period (percentile 95) | 0.10 s time step dt | 0.0045 s location (z coordinate) | 2.03 km number of elements | 10640 number of points | 10867 cost factor (nelem / dt) | 2.39e+06 max edge aspect ratio | 2.24 max equiangular skewness | 0.54 ============================================================================== GLOBAL VARIABLES: dt | 0.00446 minimum_period | 0.09995 ============================================================================== ELEMENTAL FIELDS: RHO_0 RHO_1 RHO_2 RHO_3 VP_0 VP_1 VP_2 VP_3 VS_0 VS_1 VS_2 VS_3 dt edge_aspect_ratio element_type equiangular_skewness fluid minimum_period ============================================================================== NODAL FIELDS: ============================================================================== SIDE SETS: x0 x1 y0 y1 ============================================================================== SUCCESSFULLY GENERATED MESH IN 0.157441 SECONDS. SAVED TO 'coupled_layered.e' (235.9 KB).  # View in paraview. !$PARAVIEW_BIN coupled_layered.e

# Running the elastic simulation

# This is a very short simulation for speed reasons. If you can spare a
# minute or two, make the end-time larger (1.0 seconds or more).
!\$SALVUS_BIN --mesh-file coupled_layered.e --model-file coupled_layered.e \
--polynomial-order 4 --dimension 2 --end-time 0.02 --cmd-line-source \
--source-temporal-type ricker --source-spatial-type scalar \
--source-location-x 6000 --source-location-y 4000 \
--source-scale 1e10,1e10,0 --source-center-frequency 5.0 \
--absorbing-boundaries x0,x1,y0 \
--save-movie --movie-file-name movie.h5 --movie-fields u_ELASTIC \
--save-frame-every 100 --save-fields u_ELASTIC,u_ACOUSTIC \
--io-sampling-rate 50

====================================================
Salvus version 0.0.1-686-g644f5f44-dirty
Floating point size: 64
Compiled for GLL orders: 4
Running on 1 rank(s).
====================================================

--save-wavefield-file not specified. Using default: wavefield.h5
--source-time-delay not set. Will auto-calculate.
Initializing problem.
Start time set to -0.311879.
Time step set to: 0.000384818 s.
Start time set to -0.311879.
Begin time loop.
Time loop progress [100%].
Time loop completed in 32.7965 seconds.
Begin post processing.
Problem complete.


Check out what the movie looks like, and feel free to play around a bit more. With the API of the mesher you can make models with topography, and curved fluid solid boundaries, but the workflow remains the same.