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"]
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
# 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");
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.