Benchmark problems in 2D and 3D for the elastic wave equation.
Branch: master
Clone or download
Fetching latest commit…
Cannot retrieve the latest commit at this time.
Permalink
Type Name Latest commit message Commit time
Failed to load latest commit information.
figs
literature
matrix_io
.gitignore
LICENSE
README.md
download_marmousi2.sh
elast_wedge.py
marmousi2.py

README.md

DOI

Two benchmark problems for the time-harmonic elastic wave equation in 2D and 3D

We present two benchmark problems for the time-harmonic elastic wave equation in 2D and 3D. The aim of this repository is to make the test cases considered in Section 5 of Baumann et al., 2017 publicly available.

Mathematical formulation

We consider the elastic wave equation in a frequency-domain formulation,

elastic wave eqn

where the unknown u is the displacement vector at the k-th frequency. The code we present allows an inhomogeneous medium, and implements a stress-free material-air boundary condition in the north, and first-order Sommerfeld radiation boundary conditions elsewhere.

A finite-element discretization is described in all detail in Section 2 of Baumann et al., 2016. The discrete problem yields,

Alt text

In contrast to the work in our publication, the resulting linear systems are here solved sequentially (over k) using python's sparse solver scipy.sparse.linalg.spsolve without preconditioning. For large problems, we recommend to set the flag storing=True which stores the discretization matrices K,C,M in a folder /matrix_io and does not solve the resulting linear system. The matrices are then stored in matrix market format, with a MATLAB interface matrix_io/matlab_io.m provided.

Alt text

The above spy plots can be obtained by setting spy=True. The left plot resembles a problem of ndims=2 and finite element degree=2. The right plot belongs to a problems of ndims=3and degree=1.

A sample of numerical test cases

2D elastic wedge problem

We define a new elastic wedge problem consisitng of three layers with varying physical parameters given in a computational domain [0,600]x[0,1000] meters. The parameters are given in the following table, and Lame parameters are computed accordingly.

parameter layer #1 layer #2 layer #3
rho [kg/m^3] 1800 2100 1950
c_p [m/s] 2000 3000 2300
c_s [m/s] 800 1600 1100

The python script elast_wedge.py with the following parameters,

python3 elast_wedge.py --ndims=2 --freq=[16.0] --dx=2.5 --dz=2.5 --plots=True --nprocs=4

yields the numerical results:

Alt text

3D elastic wedge problem

We extend the 2D wedge problem to 3D by expending all parameters in y-direction. A 3D discretization can be obtained by setting the flag ndims=3 and by specifying dy.

python3 elast_wedge.py --ndims=3 --freq=[4.0] --dx=10.0 --dy=10.0 --dz=10.0 --storing=True --nprocs=4

Alt text

Marmousi-II problem

In the folder /data the data set of the Marmousi-II problem is stored in segy format. We consider the following subset of the problem:

Alt text

Solving the Marmousi-II at freq=6 Hertz can be done via:

python3 marmousi2.py --freq=[6.0] --n_course=4 --nprocs=4

Here, we use every fourth grid point of the original problem in each spatial direction which yields dx=dz=5.0. The source is located at (Lx/3,0).

Alt text

Usage and installation

This code is purely written in python3 using standard numerical libraries such as NumPy and SciPy. For the finite element discretization we use nutils.

The following installation steps are necessary:

  • Clone this repository: git clone https://github.com/ManuelMBaumann/elastic_benchmarks.git
  • Install nutils via pip install git+https://github.com/joostvanzwieten/nutils@955bc67d219496e26b037f47855709a222850d7c
  • Run your first, low-frequency, 2D test case python3 elast_wedge.py --ndims=2 --dx=10.0 --dz=10.0 --freq=[4.0] --nprocs=4, and view your results with firefox ~/public_html/elast_wedge.py/latest/log.html.

For the Marmousi-II problem, two additional steps are required:

  • Download the Marmousi-II data set using the provided script download_marmousi2.sh [~ 450 Mb].
  • Install obspy via pip install obspy. In particular, we use the segy reader obspy.io.segy.core._read_segy for loading the data sets. Obspy requires libxml2 and libxslt.

A more detailed description on the installation of nutils can be found in this document.

Note: All plots will be saved as .png in a folder ~/public_html/elast_wedge.py/latest/. In the same folder, a file log.html contains information about the program evaluation and embeds all figures.

Declaration

The author is a PhD student in Numerical Analysis at TU Delft. My research is focused on linear solvers and preconditioning. I highly encourage experts in geophysics to comment on the numerical results and to get in touch.

References