An example with the 2015 M7.8 Nepal earthquake
Here I will walk you step by step through running the Nepal earthquake. This is one of the pre-canned examples which you can find in the MUD/examples folder. Note this is using only the high-rate GPS stations. If you want to reproduce the full results that we calculated in the Galetzka et al., 2015 paper in Science, that's more complicated because there we used InSAR. That's for another time. The model you will obtain here is different (this is inversion after all!) from the one in the paper but this is a useful end-to-end exercise.
Parameter files used for inversion have a .inv extension. This is not required I just do it so I have an easy way of identifying what is what. I have a folder separate from everything else where I keep my parameter (.inv) files. So create a folder where you will save all your parameter files and then copy the file nepal_version_2.inv.py from $MUD/examples/nepal_inverse\ into that folder. In my case all my inversion files are in /Users/dmelgar/Slip_inv/run/ this is what that folder looks like for me after copying the relevant file:
Open the nepal_version_2.inv.py file and change the variables that define where the inversion results and auxiliary files will be. This is a different folder than where you save your parameter files, like so:
home='/Users/dmelgar/Slip_inv/'
project_name='Nepal_example'
run_name='example'
home is where all the different inversions will live. project_name will be the folder name for this inversion project, and run_name will be the name of the inversion output files. After this is done set the flags to only initalize the folder structure by doing:
init=1 #Initalize project
make_green=0 #Compute GFs
make_synthetics=0 #Compute synthetics for a given model at given stations
G_from_file=0# =0 read GFs and create a new G, =1 load G from file
invert=0 # =1 runs inversion, =0 does nothing
Nothing else matters for now. Next step is to run it. Go to the terminal and run the script by typing:
python /Users/dmelgar/code/MudPy/examples/nepal_inverse/nepal_version2.inv.py
You should now see the entire inversion project folder structure. On my system it looks like this:
You need 3 auxiliary files to run the inversion (i) a fault geometry file (ii) a layered Earth structure file and (iii) a Green's function list file (GFlist) that specifies the station locations and types of data. MudPy requires the auxiliary files to be in specific locations so from $MUD/examples/nepal_inverse/ copy the following files:
(i) nepal_10.fault is the fault geometry and should be placed in home/project_name/data/model_info/
(ii) avouac.mod is the Earth structure model and should be placed in home/project_name/structure/
(iii) nepal_simple_v2.gflist is the GFlist file and should go in home/project_name/data/station_info/
Now you need data to invert! I've provided 5Hz GPS waveforms in $MUD/examples/nepal_inverse/data/. There is data for 5 sites and three components of motion. These are provided as SAC files. Copy all the data into home/project_name/data/waveforms
You need to edit the file by hand. The GFlist file should reflect the type of data that each waveform corresponds to (displacement, velocity, tsunami, etc.) and where the data is. The GFlist file for my system looks like:
#station lon lat "static,disp,vel,tsun,insar" "Static file, displacement file, velocity file,tsunami file,insar file" "static sigmas(n,e,u),displacement sigmas(n,e,u),velocity sigmas(n,e,u),tsunami sigma,tsunami sigma"
KKN4 85.2788 27.8008 0 1 0 0 0 /foo/bar /Users/dmelgar/Slip_inv/Nepal_example/data/waveforms/KKN4.disp /foo/bar /foo/bar /foo/bar 1 1 1 1.0 1.0 1.0 1 1 1 0 0 0 0 1 0
NAST 85.3277 27.6567 0 1 0 0 0 /foo/bar /Users/dmelgar/Slip_inv/Nepal_example/data/waveforms/NAST.disp /foo/bar /foo/bar /foo/bar 1 1 1 1.0 1.0 1.0 1 1 1 0 0 0 0 1 0
RMTE 86.5971 26.991 0 1 0 0 0 /foo/bar /Users/dmelgar/Slip_inv/Nepal_example/data/waveforms/RMTE.disp /foo/bar /foo/bar /foo/bar 1 1 1 1.0 1.0 1.0 1 1 1 0 0 0 0 1 0
SNDL 85.7989 27.3849 0 1 0 0 0 /foo/bar /Users/dmelgar/Slip_inv/Nepal_example/data/waveforms/SNDL.disp /foo/bar /foo/bar /foo/bar 1 1 1 1.0 1.0 1.0 1 1 1 0 0 0 0 1 0
SYBC 86.7124 27.8143 0 1 0 0 0 /foo/bar /Users/dmelgar/Slip_inv/Nepal_example/data/waveforms/SYBC.disp /foo/bar /foo/bar /foo/bar 1 1 1 1.0 1.0 1.0 1 1 1 0 0 0 0 1 0
CHLM 85.3141 28.2072 0 1 0 0 0 /foo/bar /Users/dmelgar/Slip_inv/Nepal_example/data/waveforms/CHLM.disp /foo/bar /foo/bar /foo/bar 1 1 1 1.0 1.0 1.0 1 1 1 0 0 0 0 1 0
The five numbers 0 1 0 0 0 indicate the data type. These stations are all labeled as displacement waveforms. One station can have multiple data types. There are 5 columns that specify where the waveform (or static data) for that station lives. You'll note 4 of those columns are just empty place holders labeled /foo/bar/ this is just to indicate there's no data for that data type. Change the displacement column (as above) to point to where you are keeping the waveform data.
The parameter file should work as described. It will calculate 5Hz Green functions and synthetics and lowpass filter them to 0.5Hz. It will use Laplacian spatial smoothing (regularization) and it will try 20 values of regularization strenght (the regularization parameter) logarithmically spaced between 10^-2 and 10^1. The epicenter are already defined and it will assume a rupture velocity of 3.2km/s and use only a single time window. These are all aprameters to be changed and explored to find better fits to the data but for now do not change them.
In preparation to run the inversion the path of least resistance is to set ncpus=1 this will use only 1 core on your machine and make the ivnersion a bit slow but should work fine. If you get tired of the inversion being slow take a look at this wiki entry on how to enable parallel processing. Most modern laptops, for example, have 4-12 cores that you can use all at once.