Mondaic - Full Waveform Solutions

This tutorial is presented as Python code running inside a Jupyter Notebook, the recommended way to use Salvus. To run it yourself you can copy/type each individual cell or directly download the full notebook, including all required files.

Download Jupyter Notebook [133.66 kB]

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.

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.

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.

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.

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

Copy

import os
import pyasdf

from salvus.flow import api, simple_config
from salvus.mesh import simple_mesh

SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")

Copy

# 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.min_period_in_seconds = 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.basic.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.730000e22,
    mtt=-2.810000e21,
    mpp=-1.450000e22,
    mrt=2.120000e22,
    mrp=4.550000e22,
    mtp=-6.570000e21,
    source_time_function=simple_config.stf.GaussianRate(
        half_duration_in_seconds=period / 2
    ),
)

# 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)

# Here we use a pre-compiled STATIONXML.
# Alternatively, we could use obspy's FDSN client to directly
# 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"
# )
w.add_receivers(
    simple_config.receiver.seismology.parse(
        "data/inventory.xml", dimensions=3, fields=["displacement"]
    )
)

# Visualize it.
w

<salvus.flow.simple_config.simulation.waveform.Waveform object at 0x7fc9af5c5310>

Copy

# 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=SALVUS_FLOW_SITE_NAME,
    output_folder="global_simulation",
)

SalvusJob `job_2403150935814350_a2e6a34353` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 0.12.16
  * Floating point size: 32

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

<salvus.flow.sites.salvus_job.SalvusJob at 0x7fc9a9086110>

Copy

# Now we'll just randomly select a waveform to plot.
ds = pyasdf.ASDFDataSet("./global_simulation/receivers.h5")
ds.waveforms.IU_ANMO.displacement.plot()