import os
import salvus_seismo
# Paths are grabbed from environment variables.
PARAVIEW_BIN = os.environ["PARAVIEW_BIN"]
SALVUS_BIN = os.environ["SALVUS_BIN"]
In this tutorial we'll run a sesimic simulation on a 2D global Earth model. We'll also use the salvus_seismo
subpackage to make adding source and receiers a breeze. While we're using a 2D example here so it can be run locally, the steps are the same for running 3D global simulations. First, let's generate a mesh for an Earth model suitable for a 200 second simulation. Use salvus_mesher
's command line interface to accomplish this.
!python -m pymesher.interface Circular2D --basic.period=200 \
--chunk2D.max_colatitude=180 --basic.model prem_iso_one_crust --overwrite
Setting up background model and element sizes... Creating the skeleton... Creating the unstructured mesh... attaching elastic parameters attaching attenuation parameters Computing mesh quality... Writing mesh to file... ============================================================================== SUMMARY OF MESH PROPERTIES: model name | prem_iso_one_crust dominant period input | 200.00 s elements per wavelength | 2.00 Courant Number | 0.60 resolved period (global max) | 195.44 s location (z coordinate) | 6356.88 km resolved period (percentile 95) | 195.44 s time step dt | 2.5241 s location (z coordinate) | 6356.88 km number of elements | 1632 number of points | 1697 cost factor (nelem / dt) | 6.47e+02 max edge aspect ratio | 12.82 max equiangular skewness | 0.62 ============================================================================== GLOBAL VARIABLES: dt | 2.52414 f_max | 1.00000 f_min | 0.00100 f_ref | 1.00000 minimum_period | 195.44021 nr_lin_solids | 5.00000 radius | 6371000.00000 w_0 | 0.00917 w_1 | 0.07173 w_2 | 0.37598 w_3 | 1.66061 w_4 | 8.66407 y_0 | 1.69035 y_1 | 1.13577 y_2 | 0.98728 y_3 | 0.92424 y_4 | 1.46135 ============================================================================== ELEMENTAL FIELDS: QKAPPA_0 QKAPPA_1 QKAPPA_2 QKAPPA_3 QMU_0 QMU_1 QMU_2 QMU_3 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 external fluid minimum_period ============================================================================== NODAL FIELDS: ============================================================================== SIDE SETS: r1 ============================================================================== SUCCESSFULLY GENERATED MESH IN 0.293197 SECONDS. SAVED TO 'Circular2D_prem_iso_one_crust_200.e' (189.6 KB).
Now, let's visualize the mesh in Paraview and inspect it.
!$PARAVIEW_BIN ./Circular2D_prem_iso_one_crust_200.e
Below, we'll rely on salvus_seismo
to generate a call to Salvus for us. This is a python
package which makes it simple to generate sources and receivers for use in global 2D/3D Earth models. No more messing around on the command line line as before! See an example below for putting an explosive source near the surface. Check out the comments below for an explanation of the individual parameters.
import salvus_seismo
# Add a source 1 meter under the surface. Specialize a 2-D "moment tensor", with
# a dominant period of 250 seconds.
src = salvus_seismo.Source(longitude=0.0, depth_in_m=1,
m_rr=1e20, m_rp=0, m_pp=1e20,
center_frequency=0.004)
# Place a single receiver at a lonitude of 45 degrees. In the 2-D globe, the
# latitude component is zero. Give it a made-up station name.
rec = salvus_seismo.Receiver(longitude=45, depth_in_m=0,
network="XX", station="AA")
# Generate the configuration object for salvus_seismo. Note the presence of "salvus_call"
# here. This can be use to pass any special commands, such as mpirun.
# This is a short run for speed reasons - choose a larger endtime for nicer
# visualizations!
config = salvus_seismo.Config(
mesh_file="./Circular2D_prem_iso_one_crust_200.e",
end_time=200,
salvus_call=SALVUS_BIN,
movie_file_name="test_movie.h5",
movie_fields=["u_ELASTIC"],
polynomial_order=4,
dimensions=2,
verbose=False)
# Ensure a clean directory. Salvus_seismo will fail to produce the configuration files if
# the directory already exists.
!rm -rf basic_example/
# Automatically generate the command line call for salvus. This takes the mesh file as an
# argument, as salvus_seismo will ensure that the source / receiver positions are adjusted
# to respect the mesh topology.
salvus_seismo.generate_cli_call(
source=src, receivers=[rec], config=config,
output_folder="basic_example",
exodus_file="./Circular2D_prem_iso_one_crust_200.e")
Wrote 1 receivers into the TOML file.
This call should generate a directory with some parameter files inside. All the receiver information is now contained within the file receivers.toml
, and a snippet containing the command line parameters can be seen in run_salvus.sh
. We'll take a look at them below.
!echo "\n\nreceivers.toml:\n"
!cat basic_example/receivers.toml
!echo "\n\nSalvus run command:\n"
!cat basic_example/run_salvus.sh
receivers.toml: [[receiver]] network = "XX" station = "AA" location = "" physical_latitude = 0.00000 physical_longitude = 45.00000 physical_depth_in_meters = 0.00000 medium = "solid" salvus_coordinates = [4504977.30294, 4504977.30294] transform_matrix = [ [0.70711, 0.70711], [-0.70711, 0.70711] ] Salvus run command: /home/boehm/SalvusInc/salvus_wave/build/salvus --dimension 2 --mesh-file ./Circular2D_prem_iso_one_crust_200.e --model-file ./Circular2D_prem_iso_one_crust_200.e --end-time 200.0 --polynomial-order 4 --receiver-file-name receiver.h5 --receiver-fields u_ELASTIC --save-movie --movie-file-name test_movie.h5 --movie-fields u_ELASTIC --receiver-toml basic_example/receivers.toml --source-toml basic_example/source.toml
This run command should contain everything you need to run Salvus
. You can go ahead and run it as below. Unfortunately, the progress reports in the notebook doesn't seem to work properly (the output is not captured until after the code is finishsed running). You can see the output if you run sh ./basic_examples/run_salvus
from the command line. Either way, this will run the 2D global simulation.
!sh ./basic_example/run_salvus.sh
==================================================== Salvus version 0.0.1-974-g7237b8b Floating point size: 64 Compiled for GLL orders: 4 Running on 1 rank(s). ==================================================== Initializing problem. Start time set to -389.848401. Time step set to: 0.217925 s. Start time set to -389.848401. Begin time loop. Time loop progress [100%]. Time loop completed in 13.8908 seconds. Begin post processing. Problem complete.
When the simulation is done, we can go ahead and view the movie using the topology output by petsc.
!python ./petsc_gen_xdmf.py ./test_movie.h5
!$PARAVIEW_BIN ./test_movie.xmf
You don't have to stop here. The source and receivers can parse any ObsPy objects and thus enable a data-driven generation of the Salvus input files. Aside from being convenient this also avoid common pitfalls like wrong unit conversions and manual human intervention.
This example directly acquires stations coordinates from the IRIS data center. Note that in 2D the latitude will be ignored and the simulation thus happens on the equatorial plane. Moment tensors and vectorial sources will consequently ignore the t (theta, south) components.
In this case the generate_cli_call()
takes an additional exodus_file
argument that takes a mesh. If given, the receivers will be placed exactly relative to the actual mesh surface.
Also note that usage in 3D is completely identical.
from obspy.clients.fdsn import Client
src = salvus_seismo.Source(longitude=0.0, depth_in_m=1,
m_rr=1e20, m_rp=0, m_pp=1e20,
center_frequency=0.004)
c = Client("IRIS")
recs = salvus_seismo.Receiver.parse(c.get_stations(network="IU"))
# Again a short simulation for speed reasons.
config = salvus_seismo.Config(
mesh_file="./Circular2D_prem_iso_one_crust_200.e",
end_time=200,
salvus_call=SALVUS_BIN,
polynomial_order=4,
verbose=False,
dimensions=2)
# Ensure a clean directory.
!rm -rf webservice_example/
salvus_seismo.generate_cli_call(
source=src, receivers=recs, config=config,
output_folder="webservice_example",
exodus_file="./Circular2D_prem_iso_one_crust_200.e")
Wrote 93 receivers into the TOML file.
!echo "\n\nreceivers.toml:\n"
!head -n 10 webservice_example/receivers.toml
!echo "\n\nSalvus run command:\n"
!cat webservice_example/run_salvus.sh
receivers.toml: [[receiver]] network = "IU" station = "AAE" location = "" physical_latitude = 0.00000 physical_longitude = 38.76560 physical_depth_in_meters = 0.00000 medium = "solid" salvus_coordinates = [4966542.42015, 3988295.19957] transform_matrix = [ [0.77971, 0.62614], Salvus run command: /home/boehm/SalvusInc/salvus_wave/build/salvus --dimension 2 --mesh-file ./Circular2D_prem_iso_one_crust_200.e --model-file ./Circular2D_prem_iso_one_crust_200.e --end-time 200.0 --polynomial-order 4 --receiver-file-name receiver.h5 --receiver-fields u_ELASTIC --receiver-toml webservice_example/receivers.toml --source-toml webservice_example/source.toml
# You can uncomment this line to visualize the receivers on the mesh in Paraview.
# !$PARAVIEW_BIN ./webservice_example/receivers_paraview.csv
!sh ./webservice_example/run_salvus.sh
==================================================== Salvus version 0.0.1-974-g7237b8b Floating point size: 64 Compiled for GLL orders: 4 Running on 1 rank(s). ==================================================== Initializing problem. Start time set to -389.848401. Time step set to: 0.217925 s. Start time set to -389.848401. Begin time loop. Time loop progress [100%]. Time loop completed in 12.4599 seconds. Begin post processing. Problem complete.
This has now generated an hdf5
file which can be opened and visualized with the ASDF sextant. Use it to open receivers.h5
.