Skip to content

Global Wave Propagation

Jupyter Notebook

This page is automatically generated from a 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.

Global Seismic Wave Propagation

This short tutorial demonstrates how to perform a global seismic wave simulation on Earth and how to work with seismological community data formats within Salvus.

A note of caution

Fully 3-D global seismic waveform simulations can be amongst the biggest numerical simulations. Across all sciences. In many cases you might not have enough computational power available to run things you would like to run. Seismologists spent many decades on developing approximate methods using various simplifications and it might be worthwhile searching the literature for these.

That being said - there is of course merit in these simulations and performing these is one of the reasons why Salvus was originally created. Salvus itself supports a number of these but this is not the topic of this tutorial.

How big is too big?

A good rule of thumb is to try to use about 5000 4th order elements in 3-D per rank/process. Thus a simulation with 1 million elements already requires a machine/cluster with about 200 cores.

Another thing to note is that these simulations scale with highest simulated frequency to the 4th power. Thus doubling the frequency results in 16 times the simulation costs.

What will we be doing in this tutorial?

We will run a simulation of the 2011 Tohoku-Oki earthquake and record it at receivers from the global seismological network (GSN). Because we want to run this on a laptop we will run this at periods of 4000 seconds. This is of course completely impractical for any real world simulation but it demonstrates how to use Salvus and the computation actually finishes in a few seconds.

One of Salvus' best features is its scalability. Changing a single variable (period here) and re-running the whole notebook could also be used to run realistic and large-scale simulations.

Steps we will perform:

  • Build a cubed-sphere mesh of the whole Earth.
  • Build a source object for the Tohoku-Oki earthquake.
  • Download receiver information for GSN stations from the IRIS DMC data center.
  • Assemble all that information into a waveform simulation object.
  • Actually run that simulation.
  • Look at the data.

Additional requirements

Aside from Salvus, this tutorial requires the pyasdf and obspy Python libraries which can be installed either via pip or via conda.

%matplotlib inline
%config Completer.use_jedi = False

import obspy.clients.fdsn
import pyasdf

from salvus_flow import api, simple_config
from salvus_mesh import simple_mesh

Setup the simulation

# Controls the dominant period of the mesh and the width
# of the source time function. It is given in seconds.
period = 4000.0

# We'll first build a mesh using the simple_mesh interface.
m = simple_mesh.Globe3D()
m.basic.period = period
# At these period we don't require a crust. Adding a 3D model
# is the topic of another tutorial.
m.basic.model = "prem_ani_no_crust"
# Higher order shapes and models better approximate the sphere.
# With order 4 we achieve a very good approximation of it
# even with only very few elements.
m.advanced.tensor_order = 4
# In order to make it a bit more interesting we'll create an
# elliptic mesh. This is the WGS84 ellipsoid.
m.spherical.ellipticity = 0.0033528106647474805
# This is an important setting. The more elementes per wavelength
# the more accurate the solution. 2 is a conservative value and
# the default. Many global seismologist only use 1 element per
# wavelength which ends up being 16 times cheaper in terms of
# simulation cost but is still usable in many scenarios.
m.advanced.elements_per_wavelength = 2.0


# Tohuko-Oki earthquake. Information is taken from the GCMT catalog
# which unfortunately does not offer a proper web service.
source = simple_config.source.seismology.SideSetMomentTensorPoint3D(
    latitude=37.5200,
    longitude=143.0500,
    depth_in_m=20000,
    side_set_name="r1",
    mrr=1.730000e+22,
    mtt=-2.810000e+21,
    mpp=-1.450000e+22,
    mrt=2.120000e+22,
    mrp=4.550000e+22,
    mtp=-6.570000e+21,
    source_time_function=simple_config.stf.GaussianRate(
        half_duration_in_seconds=period / 2)
)
# Download GSN stations via IRIS. _GSN is the virtual GSN network
# which groups a number of actual seismic network.
inv = obspy.clients.fdsn.Client("IRIS").get_stations(
    network="_GSN", level="station", format="text")


# Create the simulation object and combine all the information.
w = simple_config.simulation.Waveform(mesh=m.create_mesh())
# Sources and receivers will be placed exactly relative to the
# local mesh surface. Please refer to the sources and receivers
# documentation for more details.
w.add_sources(source)
w.add_receivers(
    simple_config.receiver.seismology.parse(
        inv, dimensions=3, fields=["displacement"]))

# Visualize it.
w

/miniconda/envs/py37/lib/python3.7/site-packages/pyasdf/asdf_data_set.py:52: UserWarning: If you install the `joblib` module, the src/rec location process would run in parallel.
  closure_warn(self, *args, **kwargs)

<salvus_flow.simple_config.simulation.Waveform object at 0x7fa8a67cf390>

Run the simulation and visualize the waveforms

# We use SalvusFlow to run the simulation. The site determines
# where it will run in the end. Might be the local machine, or
# a large remote cluster.
api.run(input_file=w, site_name="local",
        output_folder="global_simulation")
Job `job_1908041558981269_e516ddb7cb` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 0.9.8-32-g18bc6850
  * Floating point size: 32

* Downloaded 35.1 MB of results to `global_simulation`.
* Total run time: 17.05 seconds.
* Pure simulation time: 15.38 seconds.





<salvus_flow.sites.salvus_job.SalvusJob at 0x7fa8865fe320>
# Now we'll just randomly select a waveform to plot.
ds = pyasdf.ASDFDataSet("./global_simulation/receivers.h5")
ds.waveforms.IU_ANMO.displacement.plot();

png