Global 2D Simulations

Download Jupyter Notebook

This page is automatically generated from a Jupyter Notebook. You can view and download it here and here.

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.

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...
Computing resolved period...
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
  fluid
  minimum_period
==============================================================================
NODAL FIELDS:
==============================================================================
SIDE SETS:
  r1
==============================================================================

SUCCESSFULLY GENERATED MESH IN 0.340235 SECONDS.
SAVED TO 'Circular2D_prem_iso_one_crust_200.e' (186.7 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:

/Users/lion/workspace/code/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 --cmd-line-source true --source-temporal-type ricker --source-center-frequency 0.004 --save-movie --movie-file-name test_movie.h5 --movie-fields u_ELASTIC --source-location-x 6370999.0 --source-location-y 3.90111176636e-10 --source-spatial-type moment_tensor --source-scale 1e20,1e20,0 --receiver-toml basic_example/receivers.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-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 -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 14.076 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 95 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 = "AAAA"
location = ""
physical_latitude = 0.000
physical_longitude = -176.68
physical_depth_in_meters = 0.000
medium = "solid"
salvus_coordinates = [-6359208.39, -368429.13]
transform_matrix = [ [-0.998, -0.057],


Salvus run command:

/Users/lion/workspace/code/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 --cmd-line-source true --source-temporal-type ricker --source-center-frequency 0.004 --source-location-x 6370999.0 --source-location-y 3.90111176636e-10 --source-spatial-type moment_tensor --source-scale 1e20,1e20,0 --receiver-toml webservice_example/receivers.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-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 -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.2347 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.