search

seismology / mc_kernel

Calculate seismic sensitivity kernels using Axisem
GNU General Public License v3.0
24 stars 8 forks source link

readme

MC Kernel

DOI Build Status Coverage Status

Authors:

Simon Stähler, Martin van Driel, Ludwig Auer, Kasra Hosseini, Tarje Nissen-Meyer

License:

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, version 3 of the License.

Acknowledgements:

If you use this program for your work, please cite:

Simon C. Stähler, Martin van Driel, Ludwig Auer, Kasra Hosseini, Karin Sigloch, and Tarje Nissen-Meyer: MC Kernel: Broadband Waveform Sensitivity Kernels for Seismic Tomography, Geophysical Research Abstracts, Vol. 18, EGU2016-7020-2, 2016

Download

To download the current developer version, use git clone

git clone https://github.com/seismology/mc_kernel

Prerequisite

MC Kernel needs a recent Fortran compiler (gfortran 4.8 or later), an MPI installation, the NetCDF library for accessing the wavefield files, the FFTW library for Fourier transforms. The installation of these libraries can be done by hand on desktop machines or using modules on HPC environments. The recommended submit script uses Python and needs NetCDF4

Ubuntu/Debian Linux

Since Ubuntu 14.04LTS, the system libraries can be used:

- sudo apt-get install gfortran libnetcdff5 libnetcdf-dev libfftw3-dev openmpi-bin libopenmpi-dev

MacOS X

We recommend using Homebrew to download and compile the necessary libraries. Be careful to add the --with-fortran argument while installing any library. Otherwise, only C libraries are installed. If you did install the libraries before without the --with-fortran argument, it may be necessary to remove and reinstall them.

HPC environments

Be sure to load modules for Fortran, NetCDF, FFTW. On SuperMUC, the necessary commands are

module load fortran
module load netcdf/mpi
module load fftw

AxiSEM wavefields

Using the code requires computation of a global seismic wavefield using AxiSEM. Please refer to the AxiSEM documentation on how to do this. A set of example wavefields with a dominant period of 40s can be downloaded:

wget https://www.geophysik.uni-muenchen.de/~staehler/kerner_wavefields.tar.bz2
tar -xvf kerner_wavefields.tar.bz2

Installation

Change into the download directory and copy the included template files into the main directory

cd kerner
./copy_templates.sh

The file make_mc_kerner.macros allows you to modify the compiler name and compiler flags according to your system. To compile, use

make

and afterwards 

make check

to run a set of tests (wavefields downloaded in the previous step are necessary for the tests to complete).

Usage

To run the code, a convenient submit script is available:

python submit.py -n NSLAVES JOB_NAME

This creates a run directory JOB_NAME, puts the default input files there and starts an MPI job with one master and NSLAVES worker tasks. In this simple form, it uses default values stored in the files CMTSOLUTION, receivers.dat, filters.dat and stf_20s.dat together with AxiSEM wavefields in ./wavefield/fwd and ./wavefield/bwd.
Modifying the settings for the run is possible by setting one of the options of submit.py, which can be queried by

python submit.py -h

Plot results

MC Kerner comes with a set of simple Python plotting scripts to visualize kernels and the seismograms they are based on. They can be found in the folder pyfiles and need a few extra programs to run:

To run the scripts change into the job directory

cd JOB_DIR

and run

python ../pyfiles/plot_all_seismograms.py
pvpython ../pyfiles/plot_all_kernels.py
python ../pyfiles/create_composites.py

The last script creates a composite image like the one below:

Input files

There are 3 types of inputs needed to compute kernels:

Each Receiver is described by 3 parameters:

Each kernel is described by 6 parameters

These inputs are contained in two files: receiver.dat and filter.dat

receiver.dat

Example file describing a total of 4 kernels at 2 receivers; one in 30 and one in 60 degree distance. 

2                                      # Nr of receivers
Z                                      # seismogram component
30deg  60.0000   00.0000  src1 2       # Name, lat, lon, stf name, nr of kernels
P_40s Gabor_40 CC 345.0 405.0  vp      # Kernel name, Filter name, misfit, time window begin and end, model parameter
P_30s Gabor_30 CC 345.0 390.0  vp
60deg  30.0000   00.0000  src1 2       
P_40s Gabor_40 CC 580.0 640.0  vp      
P_30s Gabor_30 CC 580.0 625.0  vp

Detailed breakdown of each line:

2                                      # Nr of receivers

How many receiver blocks will be read. If there are more receiver blocks below than this number, the remaining blocks are all ignored.

Z                                      # seismogram component

Which seismogram component do we compute the kernel on. Note that M.C. Kernel needs to be run again for a different seismogram component.

30deg  60.0000   00.0000  src1 2       # Name, lat, lon, stf name, nr of kernels
P_40s Gabor_40 CC 345.0 405.0  vp      # Kernel name, Filter name, misfit, time window begin and end, model parameter
P_30s Gabor_30 CC 345.0 390.0  vp

Receiver block. Starts with one line defining the receiver location and the number of kernels described on it. Followed by one line per kernel. Each kernel line contains the 6 entries described in the comment above. Allowed values for the misfit are:

Allowed values for the model parameter are:

filters.dat

Example filters.dat file describing 3 filters:

3
filter1 Gabor 40.0 0.5 0.0 0.0  
filter2 Butterw_HP 20.   0.0 2.0 0.0  
filter3 Butterw_BP 20.0 40.0 6.0 0.0

Allowed values for the filter type are:

All Butterworth filters can be turned into a zero-phase filter by adding _ZP to the type. Note that this doubles the filter order.

FAQs

Anaconda interference

The python package manager Anaconda installs his own version of gfortran, which interferes with the system version and the MPI libraries. A solution is to replace the compiler name FC and the MPIRUN variable in make_mc_kerner.macros with absolute paths

FC = /usr/bin/mpif90

Create new AxiSEM wavefields

Prerequisites

AxiSEM

Install the latest version of AxiSEM from github. The versions published on axisem.info are rather old.

Some more comments

Remember: To compute a kernel, two groups with a total of 6 AxiSEM runs are necessary:

AxiSEM contains a Python script to group the forward and the backward runs respectively. Note the following:

AxiSEM settings

go to your AxiSEM root directory ($AXISEM)

cd $AXISEM

and copy the template input files

./copy_templates.sh

then modify two input files.

Simulation length

In inparam_basic set the simulation length to about 100 seconds after the last kernel time window that you want to use. If you want to look at

There are a number of other settings in inparam_basic, which at this stage should be left alone.

Wavefield output settings

In inparam_advanced check these parameters 

KERNEL_WAVEFIELDS  true
KERNEL_COLAT_MIN      0.
KERNEL_COLAT_MAX    180.
KERNEL_RMIN           0.
KERNEL_RMAX        7000.

KERNEL_RMIN should be set to 0 km (or the radius of the innermost cell in your kernel mesh), KERNEL_RMAX to at least the radius of your outermost mesh cell (typically the radius of the planet). It can be bigger, so 7000 is conservative for the Earth radius in kilometers. Again: Leave the other settings untouched.

Run AxiSEM

Setting up the 6 runs is a bit involved, so there is a Python script which does most of the work for you:

python ./submit.py -h
usage: submit.py [-h] [--nrad NRAD] [--ntheta NTHETA] [-j JOB_TYPE] [-r RUN_TYPE] [--max_colat MAX_COLAT]
                 [--max_depth MAX_DEPTH] [--src_depth SRC_DEPTH] [--ncl NCL] [-w WALL_TIME] [-m MAIL_ADDRESS] [-a ACCOUNT]
                 job_name mesh_file mesh_period

Create AxiSEM run and submit job.

positional arguments:
  job_name              Job directory name. 
  mesh_file             Mesh file. Choose path to external mesh file or one
                        of the following AxiSEM internal models:
                        'prem_iso', 'prem_iso_solid', 'prem_iso_onecrust',
                        'prem_iso_light', 'prem_iso_solid_light',
                        'prem_ani', 'prem_ani_onecrust', 'prem_ani_light',
                        'ak135', 'ak135f', 'iasp91'
  mesh_period           Mesh period 

optional arguments:
  -h, --help            show this help message and exit
  --nrad NRAD           Number of radial slices (default: 1)
  --ntheta NTHETA       Number of theta slices (default: 2)
                        Set to 0 to use the maximum number possible for this
                        model and frequency
  -j JOB_TYPE, --job_type JOB_TYPE
                        Job type (local, Euler, CSCS Daint)
  -r RUN_TYPE, --run_type RUN_TYPE
                        Run type (Instaseis, MC Kernel)
                        Options: 
                           instaseis (default) 
                           mckernel_fwd
                           mckernel_bwd
  --src_depth SRC_DEPTH
                        Source depth in kilometer
  --ncl NCL             Number of coarsening layers (default: 1)
  -w WALL_TIME, --wall_time WALL_TIME
                        Wall time for the solver in hours
  -m MAIL_ADDRESS, --mail_address MAIL_ADDRESS
                        Mail adress for notifications
  -a ACCOUNT, --account ACCOUNT
                        Daint project accoun

To compute wavefields for an MC Kernel run, use the -r mckernel_fwd or -r mckernel_bwd run types. For example, to fire off runs for a Mesh with 20 second period in IASP91, that uses 8 cores (use less cores if your computer has less). 

python ./submit.py FWD_iasp91 iasp91 20 -r mckernel_fwd --ntheta 8

The first argument FWD_iasp91 is a name for the run, the second argument chooses the model iasp91, the third sets the minimum period, the fourth makes it a forward run with 4 moment sources, the last chooses 8 cores. To give you an idea: This run should take about ten minutes in total on a Macbook Pro 2019. Monitor its progress with 

tail -f ./runs/FWD_iasp91/MZZ/OUTPUT

After the MZZ run is done, check the output of MXX_P_MYY, MXY_MXX_M_MYY, MXZ_MYZ. After the 4 simulations are done, the wavefield database is automatically transformed to time-contiguous.

The backward run is triggered with

python ./submit.py BWD_iasp91 iasp91 20 -r mckernel_bwd --ntheta 8

and should take half as long.

Copy the database files

After the backward run is done, check that there are database files created

ls -lh runs/FWD_iasp91/FWD_iasp91_database/MZZ/Data

This should show a file ordered_output.nc4 with roughly 100 MB in size.

Next, copy the databases to the MC Kernel wavefield directory

cp -r runs/FWD_iasp91/FWD_iasp91_database ../mc_kernel/wavefields
cp -r runs/BWD_iasp91/BWD_iasp91_database ../mc_kernel/wavefields

You can delete the AxiSEM run dir now, if you need to save diskspace. Otherwise, you're done and can run MC Kernel!