Regional 3D 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.

This tutorial is an extension of the global 3D one to show that everything works exactly the same in 3D, this time on an example that is a bit more realistic. Additionally it showcases the ObsPy integration a bit more.

We will simulate low-frequency wave originating from the March 6th Switzerland earthquake.

%matplotlib inline

import os
import warnings
warnings.simplefilter("ignore")

import salvus_seismo

# Paths are grabbed from environment variables.
PARAVIEW_BIN = os.environ["PARAVIEW_BIN"]
SALVUS_BIN = os.environ["SALVUS_BIN"]

Get event and stations

Step 1 is to get the earthquake and the stations.

# For the event we will use the Italian event web service as they do have a
# moment tensor solution of the earthquake.
from obspy.clients.fdsn import Client
import obspy

# Initialize a client for the INGV webservice.
c = Client("INGV")

cat = c.get_events(starttime=obspy.UTCDateTime(2017, 3, 6), endtime=obspy.UTCDateTime(2017, 3, 7),
                   minmagnitude=4, maxmagnitude=5, latitude=47, longitude=8, maxradius=2)

print(cat)
cat.plot(projection="local")

event = cat[0]
1 Event(s) in Catalog:
2017-03-06T20:12:06.790000Z | +46.916,   +8.931 | 4.4 ML | manual

png

# For the stations we will use the swiss service and just get
# all stations available at the time of the event.
c = Client("ETH")
inv = c.get_stations(starttime=obspy.UTCDateTime(2017, 3, 6),
                     endtime=obspy.UTCDateTime(2017, 3, 7))
inv.plot(projection="local");

png

Create a mesh

Next step is to create a mesh.

In this tutorial we'll run a sesimic simulation on a 3D regional Earth model. We'll also use the salvus_seismo subpackage to make adding source and receiers a breeze. Let's get a suitable mesh centered around Zurich.

# This nice little tool gets you the coordinates of a place of desire.
!python -m pymesher.getcoordinates zurich
best match:  Zürich, Bezirk Zürich, Zürich, Schweiz/Suisse/Svizzera/Svizra
input for pymesher SphericalChunk3D meshes with geocentric latitude
[47.180579, 8.542322, null]
!python -m pymesher.interface SphericalChunk3D \
    --basic.period=30 \
    --basic.model prem_iso_one_crust \
    --spherical.min_radius=6200 \
    --chunk3D.max_colatitude1=2.0 \
    --chunk3D.max_colatitude2=2.0 \
    --chunk3D.euler_angles 47.176739 8.540443 0 \
    --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            |     30.00 s
  elements per wavelength          |      2.00
  Courant Number                   |      0.60

  resolved period (global max)     |     27.80 s
    location (z coordinate)        |   6358.72 km
  resolved period (percentile 95)  |     27.80 s
  time step dt                     |    2.5241 s
    location (z coordinate)        |   6358.72 km
  number of elements               |       400
  number of points                 |       605
  cost factor (nelem / dt)         |  1.58e+02

  max edge aspect ratio            |      1.82
  max equiangular skewness         |      0.00
==============================================================================
GLOBAL VARIABLES:
  dt                         |   2.52414
  ellipticity                |   0.00000
  f_max                      |   1.00000
  f_min                      |   0.00100
  f_ref                      |   1.00000
  minimum_period             |  27.79868
  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
  QKAPPA_4
  QKAPPA_5
  QKAPPA_6
  QKAPPA_7
  QMU_0
  QMU_1
  QMU_2
  QMU_3
  QMU_4
  QMU_5
  QMU_6
  QMU_7
  RHO_0
  RHO_1
  RHO_2
  RHO_3
  RHO_4
  RHO_5
  RHO_6
  RHO_7
  VP_0
  VP_1
  VP_2
  VP_3
  VP_4
  VP_5
  VP_6
  VP_7
  VS_0
  VS_1
  VS_2
  VS_3
  VS_4
  VS_5
  VS_6
  VS_7
  dt
  edge_aspect_ratio
  element_type
  equiangular_skewness
  external
  fluid
  minimum_period
==============================================================================
NODAL FIELDS:
==============================================================================
SIDE SETS:
  p0
  p1
  r0
  r1
  t0
  t1
==============================================================================

SUCCESSFULLY GENERATED MESH IN 0.376695 SECONDS.
SAVED TO 'SphericalChunk3D_prem_iso_one_crust_30.e' (264.1 KB).

Now, let's visualize the mesh in Paraview and inspect it.

!$PARAVIEW_BIN ./SphericalChunk3D_prem_iso_one_crust_30.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! Check out the comments below for an explanation of the individual parameters.

import salvus_seismo

src = salvus_seismo.Source.parse(event)
recs = salvus_seismo.Receiver.parse(inv)

# Choose a center frequency suitable for our mesh.
src.center_frequency = 1.0 / 33.0

# 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 simulation for speed reasons. Set the endtime to maybe 120 seconds
# if you have some time.
config = salvus_seismo.Config(
    mesh_file="./SphericalChunk3D_prem_iso_one_crust_30.e",
    end_time=30,
    salvus_call=SALVUS_BIN,
    movie_file_name="test_movie.h5",
    movie_fields=["u_ELASTIC"],
    polynomial_order=4,
    verbose=False,
    dimensions=3)

# Ensure a clean directory. Salvus_seismo will fail to produce the configuration files if
# the directory already exists.
!rm -rf 3d_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=recs, config=config,
    output_folder="3d_example",
    exodus_file="./SphericalChunk3D_prem_iso_one_crust_30.e")
Retrieving surface information from exodus file ... [DONE in 0.03 seconds]
Automatically determined the following absorbing boundary side sets: p0, p1, r0, t0, t1
Could not locate station 8X.CA03A on mesh. Outside of mesh? latitude =  18.42, longitude =  43.68
Could not locate station 8X.CA06A on mesh. Outside of mesh? latitude =  16.25, longitude =  42.20
Could not locate station 8X.CA01A on mesh. Outside of mesh? latitude =  19.37, longitude =  43.43
Could not locate station 8X.CA07A on mesh. Outside of mesh? latitude =  16.06, longitude =  43.40
Could not locate station 8X.CA05A on mesh. Outside of mesh? latitude =  16.18, longitude =  42.86
Could not locate station 8X.CA02A on mesh. Outside of mesh? latitude =  18.59, longitude =  42.98
Could not locate station 8X.CA09A on mesh. Outside of mesh? latitude =  17.15, longitude =  42.79
Could not locate station 8X.CA08A on mesh. Outside of mesh? latitude =  17.36, longitude =  42.58
Could not locate station 8X.CA04A on mesh. Outside of mesh? latitude =  18.10, longitude =  43.11
Could not locate station Z3.A291A on mesh. Outside of mesh? latitude =  11.88, longitude =  46.43
Could not locate station Z3.A250A on mesh. Outside of mesh? latitude =  14.91, longitude =  44.89
Could not locate station Z3.A255A on mesh. Outside of mesh? latitude =  16.05, longitude =  45.32
Could not locate station Z3.A050A on mesh. Outside of mesh? latitude =  16.53, longitude =  44.29
Could not locate station Z3.A052A on mesh. Outside of mesh? latitude =  17.51, longitude =  44.91
Could not locate station Z3.A253A on mesh. Outside of mesh? latitude =  13.82, longitude =  45.10
Could not locate station Z3.A051A on mesh. Outside of mesh? latitude =  16.91, longitude =  44.82
Could not locate station Z3.A273A on mesh. Outside of mesh? latitude =  17.82, longitude =  45.72
Could not locate station Z3.A280A on mesh. Outside of mesh? latitude =   7.91, longitude =  44.35
Could not locate station Z3.A282A on mesh. Outside of mesh? latitude =   7.61, longitude =  45.06
Could not locate station Z3.A251A on mesh. Outside of mesh? latitude =  16.27, longitude =  44.94
Could not locate station Z3.A252A on mesh. Outside of mesh? latitude =  17.85, longitude =  45.24
Could not locate station Z3.A284A on mesh. Outside of mesh? latitude =   9.38, longitude =  44.94
Could not locate station Z3.A254A on mesh. Outside of mesh? latitude =  15.80, longitude =  45.12
Could not locate station Z3.A286A on mesh. Outside of mesh? latitude =   8.83, longitude =  45.17
Could not locate station Z3.A283B on mesh. Outside of mesh? latitude =   8.29, longitude =  45.05
Could not locate station Z3.A285A on mesh. Outside of mesh? latitude =   9.91, longitude =  44.70
Could not locate station Z3.A271A on mesh. Outside of mesh? latitude =  18.83, longitude =  46.96
Could not locate station Z3.A281A on mesh. Outside of mesh? latitude =   7.71, longitude =  44.66
Could not locate station Z3.A272A on mesh. Outside of mesh? latitude =  18.97, longitude =  46.55
Could not locate station Z3.A268A on mesh. Outside of mesh? latitude =  17.92, longitude =  47.24
Could not locate station DK.KULLO on mesh. Outside of mesh? latitude = -57.22, longitude =  74.48
Could not locate station DK.ILULI on mesh. Outside of mesh? latitude = -51.10, longitude =  69.08
Could not locate station DK.NUUG on mesh. Outside of mesh? latitude = -53.20, longitude =  71.42
Wrote 183 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"
!head -n 20 3d_example/receivers.toml

!echo "\n\nSalvus run command:\n"
!cat 3d_example/run_salvus.sh
receivers.toml:

[[receiver]]
network = "C4"
station = "CERN5"
location = ""
physical_latitude = 46.11764
physical_longitude = 6.07608
physical_depth_in_meters = 0.00000
medium = "solid"
salvus_coordinates = [4391392.73248, 467450.31225, 4591941.82537]
transform_matrix = [ [0.68929, 0.07337, 0.72076],
                     [-0.71672, -0.07629, 0.69318],
                     [-0.10585, 0.99438, 0.00000] ]

[[receiver]]
network = "C4"
station = "CERNS"
location = ""
physical_latitude = 46.07468
physical_longitude = 6.06621
physical_depth_in_meters = 0.00000


Salvus run command:

/home/boehm/SalvusInc/salvus_wave/build/salvus --dimension 3 --mesh-file ./SphericalChunk3D_prem_iso_one_crust_30.e --model-file ./SphericalChunk3D_prem_iso_one_crust_30.e --end-time 30.0 --polynomial-order 4 --receiver-file-name receiver.h5 --receiver-fields u_ELASTIC --absorbing-boundaries p0,p1,r0,t0,t1 --save-movie --movie-file-name test_movie.h5 --movie-fields u_ELASTIC --receiver-toml 3d_example/receivers.toml --source-toml 3d_example/source.toml
!sh ./3d_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 -51.459989.
Time step set to: 0.281515 s.
Start time set to -51.459989.
Begin time loop.
Time loop progress [100%].
Time loop completed in 2.94579 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

Feel free to load the data in the ASDF sextant to view it in a convenient way.