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.

%matplotlib inline

# Import modules and setup paths.
import os
import h5py
import matplotlib.pyplot as plt
import pyasdf

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

MESH2D = "Quad_IsotropicElastic2D_Nodal_44x44.e"
MESH3D = "Hex_IsotropicElastic3D_Nodal_7x7x7.e"

Advanced input and output

In this example we will use advanced options to specify sources and receivers and explore various options for wavefield data output. We start with a very simple example of a single source-receiver pair in a homogeneous medium. The path to the mesh is also stored in the HOMOGENEOUS_MESH variable defined above. Open the file with Paraview and explore its contents.

!$PARAVIEW_BIN $MESH2D

You should already be familiar with Paraview, so just make sure that the model is indeed homogeneous by checking the values of RHO, LAMBDA, MU, VP and VS. Why are there 5 different parameters you ask? Good question. Of course, in an isotropic medium we have the relations

\[ v_p = \sqrt{\frac{\lambda + 2\mu}{\rho}},\quad\quad v_s = \sqrt{\frac{\mu}{\rho}}, \]

but for convenience we added the other fields as well. Let's check whether the values are consistent by using the Calculator filter in Paraview. From the top menu select Filters -> Data Analysis -> Calculator. In the panel on the left you can now define pointwise operations on the parameter fields. For instance, type VS_check in the field Result Array Name and put sqrt(MU/RHO) in the line below. After clicking on Apply you can now select VS_check from the dropdown menu on the top an plot the new field. The values should be the same as in VS. Repeat the same steps to check whether VP is also consistent with LAMBDA and MU.

Next, we need to decide where to put the source and the receiver in this example. Scroll down the property panel and tick the checkbox Axes Grid. You will see that our domain is given by [-1000,1000] x [0,2000]. So let's put the source at [-500,1000] and the receiver at [500,1000].

Source and receiver input

You can specify almost all input parameters using the command line. However, this gets really tedious when dealing with many sources and receivers. Fortunately, there is another way using a toml file. toml is a minimal configuration file format that is easy to read and to parse. A receiver needs at least to set the network and station name following the SEED format (http://www.fdsn.org/seed_manual/SEEDManual_V2.4.pdf). Furthermore, we need to specify whether the receiver is in fluid or solid media and its location, of course. This is how the toml file can be created:

with open("salvus.toml", "w") as text_file:
    recLocX = 500
    recLoxY = 1000
    print("\n[[receiver]]", file=text_file)
    print("network = \"AB\"", file=text_file)
    print("station = \"CED\"", file=text_file)
    print("medium = \"solid\"", file=text_file)
    print("salvus_coordinates = [{0}, {1}]".format(recLocX, recLoxY), file=text_file)

    recLocX = 600
    recLoxY = 1000
    print("\n[[receiver]]", file=text_file)
    print("network = \"AB\"", file=text_file)
    print("station = \"CEF\"", file=text_file)
    print("medium = \"solid\"", file=text_file)
    print("salvus_coordinates = [{0}, {1}]".format(recLocX, recLoxY), file=text_file)
! cat salvus.toml
[[receiver]]
network = "AB"
station = "CED"
medium = "solid"
salvus_coordinates = [500, 1000]

[[receiver]]
network = "AB"
station = "CEF"
medium = "solid"
salvus_coordinates = [600, 1000]

Check the contents of the file and make sure it contains the correct values.

Now it's time to run a simulation.

!$SALVUS_BIN --mesh-file $MESH2D --model-file $MESH2D \
 --polynomial-order 4 --dimension 2 --end-time 1.04 \
 --cmd-line-source --source-temporal-type ricker \
 --source-spatial-type vector \
 --source-location-x -500 --source-location-y 1000 \
 --source-scale 1e9,1e9 --source-center-frequency 14.5 \
 --receiver-toml salvus.toml --receiver-file-name receiver.h5 \
 --receiver-fields u_ELASTIC
====================================================
Salvus version 0.0.1-686-g644f5f44-dirty
Floating point size: 64
Compiled for GLL orders: 4
Running on 1 rank(s).
====================================================

--source-time-delay not set. Will auto-calculate.
Initializing problem.
Start time set to -0.107544.
Time step set to: 0.000735822 s.
Start time set to -0.107544.
Begin time loop.
Time loop progress [100%].
Time loop completed in 7.18417 seconds.
Begin post processing.
Problem complete.

Did it work? If you are not sure check the contents of receiver.h5:

! h5dump -H receiver.h5
HDF5 "receiver.h5" {
GROUP "/" {
   ATTRIBUTE "file_format" {
      DATATYPE  H5T_STRING {
         STRSIZE 4;
         STRPAD H5T_STR_NULLPAD;
         CSET H5T_CSET_ASCII;
         CTYPE H5T_C_S1;
      }
      DATASPACE  SCALAR
   }
   ATTRIBUTE "file_format_version" {
      DATATYPE  H5T_STRING {
         STRSIZE 5;
         STRPAD H5T_STR_NULLPAD;
         CSET H5T_CSET_ASCII;
         CTYPE H5T_C_S1;
      }
      DATASPACE  SCALAR
   }
   GROUP "Waveforms" {
      GROUP "AB.CED" {

         DATASET "AB.CED..FHE__1970-01-01T00:00:00__1970-01-01T00:00:01__displacement" {
            DATATYPE  H5T_IEEE_F64LE
            DATASPACE  SIMPLE { ( 1561 ) / ( 1561 ) }
            ATTRIBUTE "sampling_rate" {
               DATATYPE  H5T_IEEE_F64LE
               DATASPACE  SCALAR
            }
            ATTRIBUTE "starttime" {
               DATATYPE  H5T_STD_I64LE
               DATASPACE  SCALAR
            }
         }

         DATASET "AB.CED..FHZ__1970-01-01T00:00:00__1970-01-01T00:00:01__displacement" {
            DATATYPE  H5T_IEEE_F64LE
            DATASPACE  SIMPLE { ( 1561 ) / ( 1561 ) }
            ATTRIBUTE "sampling_rate" {
               DATATYPE  H5T_IEEE_F64LE
               DATASPACE  SCALAR
            }
            ATTRIBUTE "starttime" {
               DATATYPE  H5T_STD_I64LE
               DATASPACE  SCALAR
            }
         }
         DATASET "StationXML" {
            DATATYPE  H5T_STD_I8LE
            DATASPACE  SIMPLE { ( 715 ) / ( 715 ) }
         }
      }
      GROUP "AB.CEF" {

         DATASET "AB.CEF..FHE__1970-01-01T00:00:00__1970-01-01T00:00:01__displacement" {
            DATATYPE  H5T_IEEE_F64LE
            DATASPACE  SIMPLE { ( 1561 ) / ( 1561 ) }
            ATTRIBUTE "sampling_rate" {
               DATATYPE  H5T_IEEE_F64LE
               DATASPACE  SCALAR
            }
            ATTRIBUTE "starttime" {
               DATATYPE  H5T_STD_I64LE
               DATASPACE  SCALAR
            }
         }

         DATASET "AB.CEF..FHZ__1970-01-01T00:00:00__1970-01-01T00:00:01__displacement" {
            DATATYPE  H5T_IEEE_F64LE
            DATASPACE  SIMPLE { ( 1561 ) / ( 1561 ) }
            ATTRIBUTE "sampling_rate" {
               DATATYPE  H5T_IEEE_F64LE
               DATASPACE  SCALAR
            }
            ATTRIBUTE "starttime" {
               DATATYPE  H5T_STD_I64LE
               DATASPACE  SCALAR
            }
         }
         DATASET "StationXML" {
            DATATYPE  H5T_STD_I8LE
            DATASPACE  SIMPLE { ( 715 ) / ( 715 ) }
         }
      }
   }
}
}

Well, at least the traces are there. But better to plot the traces instead. Try it yourself using pyasdf.

with pyasdf.ASDFDataSet("receiver.h5") as ds:
    plt.plot(ds.waveforms.AB_CED.displacement[0].data)

png

Can you explain the different wave fronts you see in the seismgram? Where do they come from?

To verify your assumption we will visualize the entire wavefield in the next section. But before we continue, let's reduce the number of command line options for the source.

Similarly to the receivers, we can define one or more sources in a toml file instead of passing all the options to the command line. Salvus can transform spherical coordinates into the internally used cartesian coordinates, but this currently only works for receivers and not for sources. Thus, the toml field for sources is called just location instead of salvus_coordinates.

with open("salvus.toml", "a") as text_file:
    name = "source1"
    srcLocX = -500.0
    srcLocY = 1000.0
    spatial_type = "vector"
    temporal_type = "ricker"
    scaleX = 1e9
    scaleY = 1e9
    center_frequency = 14.5
    print("\n[[source]]", file=text_file)
    print("name = \"{0}\"".format(name), file=text_file)
    print("location = [{0}, {1}]".format(srcLocX, srcLocY), file=text_file)
    print("spatial_type = \"{0}\"".format(spatial_type), file=text_file)
    print("temporal_type = \"{0}\"".format(temporal_type), file=text_file)
    print("center_frequency = {0}".format(center_frequency), file=text_file)
    print("scale = [{0}, {1}]".format(scaleX,scaleY), file=text_file)

Note that we append the sources to the same toml file, so be careful when you run this cell several times.

!$SALVUS_BIN --mesh-file $MESH2D --model-file $MESH2D \
 --polynomial-order 4 --dimension 2 --end-time 1.04 \
 --source-toml salvus.toml \
 --receiver-toml salvus.toml --receiver-file-name receiver1.h5 \
 --receiver-fields u_ELASTIC
====================================================
Salvus version 0.0.1-686-g644f5f44-dirty
Floating point size: 64
Compiled for GLL orders: 4
Running on 1 rank(s).
====================================================

Initializing problem.
Start time set to -0.107544.
Time step set to: 0.000735822 s.
Start time set to -0.107544.
Begin time loop.
Time loop progress [100%].
Time loop completed in 7.23567 seconds.
Begin post processing.
Problem complete.

Plot and compare both receiver.

with pyasdf.ASDFDataSet("receiver.h5") as ds, \
        pyasdf.ASDFDataSet("receiver1.h5") as ds1:
 plt.plot(ds.waveforms.AB_CED.displacement.select(component="Z")[0].data, color="green")
 plt.plot(ds1.waveforms.AB_CED.displacement.select(component="Z")[0].data, color="red")

png

The traces should be exactly the same.

Wavefield output

Volume data

So far, we have generated receiver data with Salvus and visualized the wavefield using the --save-movie option. Now, it is time to learn another way to output general wavefield data with Salvus.

--save-movie has two limitations: You can only output the wavefield at the vertices of the elements, but not at all the GLL points, which means that you are missing 117 of 125 points for a 4th-order hexahedral. Second, --save-movie is restricted to output only one field at the time. Instead, we will use the option --save-fields now.

With --save-fields you can output a variety of different fields (displacement, velocities, ...) to an hdf5 file and Salvus will automatically create an xdmf file for you to visualize the data with Paraview. So no need to call petsc_gen_xdmf.py anymore. However, one drawback at this stage is that Salvus only outputs the elemental data at all GLL points, but no topology. So you need to apply an interpolation filter (i.e. Delaunay 2D/3D) to obtain a continous field.

In addition, we need to define the output file using --save-wavefield-file (e.g., wavefield.h5) and, optionally, the sampling rate with --io-sampling-rate.

!$SALVUS_BIN --mesh-file $MESH2D --model-file $MESH2D \
 --polynomial-order 4 --dimension 2 --end-time 1.04 \
 --source-toml salvus.toml \
 --receiver-toml salvus.toml \
 --receiver-file-name receiver.h5 --receiver-fields u_ELASTIC \
 --save-fields u_ELASTIC --save-wavefield-file "wavefield.h5" --io-sampling-rate 10
====================================================
Salvus version 0.0.1-686-g644f5f44-dirty
Floating point size: 64
Compiled for GLL orders: 4
Running on 1 rank(s).
====================================================

Initializing problem.
Start time set to -0.107544.
Time step set to: 0.000735822 s.
Start time set to -0.107544.
Begin time loop.
Time loop progress [100%].
Time loop completed in 6.95663 seconds.
Begin post processing.
Problem complete.

Let's check the contents of the current folder. You should see the files wavefield.h5 and wavefield_ELASTIC.xdmf

! ls -l
total 241856
-rw-r--r--  1 lion  staff      87008 Jun 28 12:02 Hex_IsotropicElastic3D_Nodal_7x7x7.e
-rw-r--r--  1 lion  staff     261024 Jun 20 12:17 Quad_IsotropicElastic2D_Nodal_44x44.e
-rw-r--r--  1 lion  staff     105575 Jun 28 12:11 Salvus_IO.ipynb
-rw-r--r--  1 lion  staff      16222 Jun 27 15:40 curldiv0.png
-rw-r--r--  1 lion  staff      56304 Jun 27 15:40 curldiv1.png
-rw-r--r--  1 lion  staff      56225 Jun 27 15:40 curldiv2.png
-rw-r--r--  1 lion  staff      59704 Jun 28 12:12 receiver.h5
-rw-r--r--  1 lion  staff      61112 Jun 28 12:12 receiver1.h5
-rw-r--r--  1 lion  staff        357 Jun 28 12:11 salvus.toml
-rw-r--r--  1 lion  staff     383628 Jun 27 15:40 surface_data.pvsm
-rw-r--r--  1 lion  staff     106191 Jun 27 15:40 surface_snapshot.png
-rw-r--r--  1 lion  staff  122364744 Jun 28 12:12 wavefield.h5
-rw-r--r--  1 lion  staff     155686 Jun 28 12:12 wavefield_ELASTIC.xdmf
-rw-r--r--  1 lion  staff      90445 Jun 28 12:09 wavefield_ELASTIC_BND.xdmf

Now let's open the xdmf file with Paraview. You can either look at the point cloud or apply the Delaunay 2D filter Filters -> Alphabetical -> Delaunay 2D.

! $PARAVIEW_BIN wavefield_ELASTIC.xdmf

Ah, we forgot to set absorbing boundaries! That's why the seismograms looked messy. Try again with absorbing boundaries everywhere and set the sampling rate to 1.

!$SALVUS_BIN --mesh-file $MESH2D --model-file $MESH2D \
 --polynomial-order 4 --dimension 2 --end-time 1.04 \
 --source-toml salvus.toml \
 --receiver-toml salvus.toml --receiver-file-name receiver.h5 --receiver-fields u_ELASTIC \
 --absorbing-boundaries x0,x1,y0,y1 \
 --save-fields u_ELASTIC --save-wavefield-file "wavefield.h5" --io-sampling-rate 1
====================================================
Salvus version 0.0.1-686-g644f5f44-dirty
Floating point size: 64
Compiled for GLL orders: 4
Running on 1 rank(s).
====================================================

Initializing problem.
Start time set to -0.107544.
Time step set to: 0.000735822 s.
Start time set to -0.107544.
Begin time loop.
Time loop progress [100%].
Time loop completed in 9.44389 seconds.
Begin post processing.
Problem complete.

Check the wavefield and the receiver again.

!$PARAVIEW_BIN wavefield_ELASTIC.xdmf
with pyasdf.ASDFDataSet("receiver.h5") as ds:
    ds.waveforms.AB_CED.displacement.plot()

png

Until now, we only looked at the displacement fields. But you can output other fields easily with --save-fields, e.g., velocities (v_ELASTIC), strains (strain) or the full gradient (grad). Let's give it a try with the gradient to compute curl and divergence.

!$SALVUS_BIN --mesh-file $MESH2D --model-file $MESH2D \
 --polynomial-order 4 --dimension 2 --end-time 1.04 \
 --source-toml salvus.toml \
 --receiver-toml salvus.toml --receiver-file-name receiver.h5 --receiver-fields u_ELASTIC \
 --absorbing-boundaries x0,x1,y0,y1 \
 --save-fields grad --save-wavefield-file "wavefield.h5" --io-sampling-rate 100
====================================================
Salvus version 0.0.1-686-g644f5f44-dirty
Floating point size: 64
Compiled for GLL orders: 4
Running on 1 rank(s).
====================================================

Initializing problem.
Start time set to -0.107544.
Time step set to: 0.000735822 s.
Start time set to -0.107544.
Begin time loop.
Time loop progress [100%].
Time loop completed in 7.40487 seconds.
Begin post processing.
Problem complete.

Let's look at the contents of the h5 file.

! h5ls -rv wavefield.h5
Opened "wavefield.h5" with sec2 driver.
/                        Group
    Location:  1:96
    Links:     1
/ELASTIC                 Group
    Location:  1:800
    Links:     1
/ELASTIC/coordinates     Dataset {1936/1936, 25/25, 2/2}
    Location:  1:1832
    Links:     1
    Storage:   774400 logical bytes, 774400 allocated bytes, 100.00% utilization
    Type:      native double
/ELASTIC/data            Dataset {17/17, 1936/1936, 4/4, 25/25}
    Attribute: DIMENSION_LABELS {4}
        Type:      variable-length null-terminated ASCII string
        Data:  "time", "element", "[ grad_11 | grad_12 | grad_21 | grad_22 ]",
               "point"
    Location:  1:776872
    Links:     1
    Storage:   26329600 logical bytes, 26329600 allocated bytes, 100.00% utilization
    Type:      native double
/ELASTIC/time            Dataset {17/17}
    Location:  1:777216
    Links:     1
    Storage:   136 logical bytes, 136 allocated bytes, 100.00% utilization
    Type:      native double

Can you find the related field names?

They are called grad_11, grad_12, grad_12 and grad_22.

Now open the file in paraview and use the Calculator filter to compute curl and divergence. In case you don't remember, they are defined as follows in 2D $$ \text{div}~u = \frac{\partial u_1}{\partial x_1} + \frac{\partial u_2}{\partial x_2}, $$ $$ \text{curl}~u = \frac{\partial u_2}{\partial x_1} - \frac{\partial u_1}{\partial x_2}. $$

!$PARAVIEW_BIN wavefield_ELASTIC.xdmf

divergence and curl, snapshot1 divergence and curl, snapshot2 divergence and curl, snapshot3

By the way, you could use the same field tags for the receiver-fields as well in case you want to measure rotations, for instance.

Boundary data

In many applications, for instance, ambient noise tomography, we are only interested in the wavefield at the surface of the Earth. It would be a waste to store the wavefield everywhere and very inefficient to use a huge number of receivers for that, right?

This is why Salvus allows you to output boundary data at arbitrary boundaries. Again, you are free to choose which fields you want to plot. How about velocities?

The relevant options are: --save-boundary-fields and --save-boundaries to specify which fields to store and where.

Of course, it is more interesting to look at boundary data in 3D. So we need a new mesh for this.

This is a simple 3D mesh with 343 elements. Put a source in the center of the domain and use a simulation time of 25 s and frequency of 0.1 Hz. You are free to choose which boundary you want to save, but we recommend to use z1.

!$SALVUS_BIN --mesh-file $MESH3D --model-file $MESH3D \
 --polynomial-order 4 --dimension 3 --end-time 25 --cmd-line-source \
 --source-temporal-type ricker --source-spatial-type vector \
 --source-location-x 50000 --source-location-y 50000 --source-location-z 50000 \
 --source-scale 0,0,1e9 --source-center-frequency 0.1 \
 --absorbing-boundaries x0,x1,y0,y1,z0 \
 --save-boundary-fields v_ELASTIC --save-boundaries z1 \
 --save-fields v_ELASTIC --io-sampling-rate 5
====================================================
Salvus version 0.0.1-686-g644f5f44-dirty
Floating point size: 64
Compiled for GLL orders: 4
Running on 1 rank(s).
====================================================

--save-wavefield-file not specified. Using default: wavefield.h5
--source-time-delay not set. Will auto-calculate.
Initializing problem.
Start time set to -15.593936.
Time step set to: 0.127591 s.
Start time set to -15.593936.
Begin time loop.
Time loop progress [100%].
Time loop completed in 2.47129 seconds.
Begin post processing.
Problem complete.

Open Paraview and then click on File -> Load State and select the file surface_data.pvsm. Make sure to adjust the paths to your local computer!

!$PARAVIEW_BIN

surface data snapshot

Alright, that's it. Enough colorful pictures. Now it's time for some math.

Fin.