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:
Aside from Salvus, this tutorial requires the
obspy Python libraries which can be installed either via
pip or via
%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
# 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.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 ), ) # 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:56: 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 0x7f94d2fdcf60>
# 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_1910230016871866_cbc29710c2` running on `local` with 4 rank(s). Site information: * Salvus version: 0.10.6 * Floating point size: 32
* Downloaded 35.3 MB of results to `global_simulation`. * Total run time: 19.65 seconds. * Pure simulation time: 18.84 seconds.
<salvus_flow.sites.salvus_job.SalvusJob at 0x7f94f818ff28>
# Now we'll just randomly select a waveform to plot. ds = pyasdf.ASDFDataSet("./global_simulation/receivers.h5") ds.waveforms.IU_ANMO.displacement.plot()
<Figure size 800x750 with 3 Axes>