Documentation for computing source localization with a dense seismic array
This documentation explains how we performed the computation of matched field processing (MFP) localization with Argentière glacier RESOLVE ZO_2018 dataset.
Our method can be applied with any other dataset, and do not hesitate to contact us if you have some questions or issues. This documentation is related to this publication:
We used the process described in this documentation also in two following papers:
Observing the subglacial hydrology network and its dynamics with a dense seismic array
Here is the gitlab repository associated to the RESOLVE project.
1. Download one day of data
The first step is to download the data, which we acquired with 98 3-components seismic sensors during a one month-long experiment with a 500 Hz sampling rate.
The RESIF ZO_2018 dataset [1] should be cited as:
[1] Roux, P., Gimbert, F., & RESIF. (2021). Dense nodal seismic array temporary experiment on Alpine Glacier of Argentière (RESIF-SISMOB) [Data set]. Réseau Sismologique et géodésique Français. https://doi.org/10.15778/RESIF.ZO2018
You can click here to access the landing page where you can download the data:
ZO_2018 landing page from RESIF portal
There are many ways to do this depending on the amount of data.
1.a. Downloading with FDSN ph5ws webservices from RESIF datacenter (< 100GB data)
Here we give an example where we download 5-min of data for the 98 active sensors (only vertical component).
> wget "https://ph5ws.resif.fr/fdsnws/dataselect/1/query?net=ZO&channel=DPZ&start=2018-05-01T00:00:00Z&end=2018-05-01T00:05:00Z&format=mseed" -O data.mseed
This is an example of the ouptut message:
--2021-11-03 10:14:59-- https://ph5ws.resif.fr/fdsnws/dataselect/1/query?net=ZO&channel=DPZ&start=2018-05-01T00:00:00Z&end=2018-05-01T00:05:00Z&format=mseed
Resolving ph5ws.resif.fr (ph5ws.resif.fr)... 129.88.201.60
Connecting to ph5ws.resif.fr (ph5ws.resif.fr)|129.88.201.60|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 15310848 (15M) [application/vnd.fdsn.mseed]
Saving to: ‘data.mseed’
data.mseed 100%[============================================================>] 14.60M --.-KB/s in 0.05s
2021-11-03 10:15:07 (266 MB/s) - ‘data.mseed’ saved [15310848/15310848]
It delivers a 15MB miniSEED file. You can check the list of data you downloaded using (for example) IRIS miniSEED inspector tool msi:
> msi -tg data.mseed
Source Start sample End sample Gap Hz Samples
ZO_AR001__DPZ 2018,121,00:00:00.001999 2018,121,00:04:59.999999 == 500 150000
ZO_AR002__DPZ 2018,121,00:00:00.000000 2018,121,00:04:59.998000 == 500 150000
ZO_AR004__DPZ 2018,121,00:00:00.000000 2018,121,00:04:59.998000 == 500 150000
ZO_AR005__DPZ 2018,121,00:00:00.000000 2018,121,00:04:59.998000 == 500 150000
ZO_AR006__DPZ 2018,121,00:00:00.001998 2018,121,00:04:59.999998 == 500 150000
ZO_AR007__DPZ 2018,121,00:00:00.000000 2018,121,00:04:59.998000 == 500 150000
.
.
.
ZO_AR097__DPZ 2018,121,00:00:00.001998 2018,121,00:04:59.999998 == 500 150000
ZO_AR098__DPZ 2018,121,00:00:00.001999 2018,121,00:04:59.999999 == 500 150000
ZO_AR099__DPZ 2018,121,00:00:00.000000 2018,121,00:04:59.998000 == 500 150000
ZO_AR100__DPZ 2018,121,00:00:00.001998 2018,121,00:04:59.999998 == 500 150000
Total: 98 trace(s) with 98 segment(s)
The ph5ws webservice delivers data in .mseed format. On RESIF datacenter, these "nodes-type" data are actually stored as IRIS Passcal PH5 format. Thus the ph5ws webservice automatically converts from PH5 to miniSEED format. It is mainly because of this on-the-fly conversion that this webservice is not able to deliver a large amount of data. RESIF documentation indicates a limit around 100GB of data, but our own tests seem to point to a limit around only 1GB.
The MFP Fortran code can only read home-preprocessed HDF5 formatted datasets (that are different from IRIS PH5 format). Section 2 explains how to do that.
Note that when selecting a given period (here 5 min) the .mseed files might be cut at the first and/or last sample (for example 00:00:00.001998). We suggest you to download a bit larger time window and then cut (with using python-obspy module for example).
Il faut expliquer que dans notre travail, nous n'avions pas à gérer ces deux points : 1) premiers et derniers echantillons parfois coupés, et 2) desynchronisation d'une ou deux microsecondes pour certaines traces, car 1) on accédait directement au miniSEED complets sur disque et 2) les données sont passées par la moulinette SISMOB qui a corrigé cette desynchronisation : cette étape a été effectué par RESIF avant qu'ils ne nous mettent à disposition les miniSEED complets sur disque. On recherche actuellement trace et référence de ce traitement "moulinette SISMOB" auprès de RESIF...
1.b. Downloading with RSYNC MODULE from RESIF datacenter (> 100GB data)
This option is preferable for large dataset.
You can find the documentation on large dataset synchronization with rsync.resif.fr. Here is an example:
> DESTINATION="/tmp" # choose your local destination directory
> RSYNC_MODULE=ZO2018-PH5-DATA
> RSYNC_SERVER="rsync.resif.fr"
> RSYNC_OPTS="-rltvh --compress-level=1"
> DRYRUN="-n --stats"
> RSYNC_USERNAME=""
# only a dry-run
> rsync $RSYNC_OPTS $DRYRUN rsync://$RSYNC_USERNAME@$RSYNC_SERVER/$RSYNC_MODULE $DESTINATION
R\#303\#251sif rsync server
documentation : https://seismology.resif.fr/rsync/
support: sismo-help@resif.fr
receiving incremental file list
./
README
master.ph5
miniPH5_00001.ph5
miniPH5_00002.ph5
[...]
miniPH5_00096.ph5
miniPH5_00097.ph5
miniPH5_00098.ph5
Number of files: 101 (reg: 100, dir: 1)
[...]
Total file size: 541.28G bytes
[...]
total size is 541.28G speedup is 181,271,560.65 (DRY RUN)
# ask for total volumetry (ls -l)
> rsync rsync://$RSYNC_USERNAME@$RSYNC_SERVER/$RSYNC_MODULE
R\#303\#251sif rsync server
documentation : https://seismology.resif.fr/rsync/
support: sismo-help@resif.fr
drwxr-sr-x 12,288 2021/11/09 11:12:25 .
-rwxr-sr-x 320 2021/11/09 11:12:25 README
-rw-r--r-- 163,487,945 2020/01/25 14:50:42 master.ph5
-rw-r--r-- 6,108,098,225 2020/01/20 14:59:39 miniPH5_00001.ph5
-rw-r--r-- 5,533,694,155 2020/01/20 15:23:55 miniPH5_00002.ph5
[...]
-rw-r--r-- 5,329,278,727 2020/01/25 10:44:56 miniPH5_00096.ph5
-rw-r--r-- 5,741,636,187 2020/01/25 12:49:23 miniPH5_00097.ph5
-rw-r--r-- 5,542,145,229 2020/01/25 14:50:42 miniPH5_00098.ph5
# transfer
> rsync $RSYNC_OPTS rsync://$RSYNC_USERNAME@$RSYNC_SERVER/$RSYNC_MODULE $DESTINATION
[...]
Using RSYNC_MODULE allows to download PH5 files from RESIF datacenter. Then the user has to convert these files into miniSEED files, using the PH5 library.
2. Preprocess the data
After downloading miniSEED files from RESIF datacentre, you could use python-obspy module for preprocessing this mseed file and converting it into a HDF5 file adapted for MFP code inputs.
Once you downloaded the data in .mseed you should apply a number of steps to convert them in .hdf5 format, which is the format of input for the MFP codes we developed. We fill gaps with zero, we properly cut daily 500Hz traces (thus 43,200,200 samples), we write each daily trace as HDF5 dataset with Float32 encoding (as RESIF initially distributed us Float32 miniSEED files, that is another difference with the new PH5 format that encodes Int32 data...) and gzip compression (lossless compression which allowed us about 200% compression rate). You can find here in experiment EXP000 the description and the associated codes of all these steps.
The resulting HDF5 files (here for the vertical component) should be like this, with Float32 encoding datasets, with or without dataset GZIP-compression:
> h5ls ZO_2018_121.h5
ZO.AR001.00.DPZ.D.2018.121 Dataset {43200000}
ZO.AR002.00.DPZ.D.2018.121 Dataset {43200000}
ZO.AR004.00.DPZ.D.2018.121 Dataset {43200000}
ZO.AR005.00.DPZ.D.2018.121 Dataset {43200000}
ZO.AR006.00.DPZ.D.2018.121 Dataset {43200000}
.
.
.
ZO.AR095.00.DPZ.D.2018.121 Dataset {43200000}
ZO.AR096.00.DPZ.D.2018.121 Dataset {43200000}
ZO.AR097.00.DPZ.D.2018.121 Dataset {43200000}
ZO.AR098.00.DPZ.D.2018.121 Dataset {43200000}
ZO.AR099.00.DPZ.D.2018.121 Dataset {43200000}
ZO.AR100.00.DPZ.D.2018.121 Dataset {43200000}
> h5ls -v ZO_2018_121.h5/ZO.AR084.00.DPZ.D.2018.121
Opened "ZO_2018_121.h5" with sec2 driver.
ZO.AR084.00.DPZ.D.2018.121 Dataset {43200000/43200000}
Location: 1:6463553696
Links: 1
Chunks: {21094} 84376 bytes
Storage: 172800000 logical bytes, 78849872 allocated bytes, 219.15% utilization
Filter-0: fletcher32-3 {}
Filter-1: deflate-1 OPT {1}
Type: native float
For MFP computation, it is not mandatory to have daily HDF5 datasets, they could be hourly, several days, ... but it is important that all HDF5 datasets have the same number of samples, and are synchronized (that's why we filled gaps with zero). Indeed, these HDF5 input files don't contain any time vector dimension, all HDF5 datasets are supposed to be along the same time-vector.
Besides, here we choose the SDS naming convention for naming HDF5 datasets, but it is possible to choose other names. What's important is to provide the same HDF5 dataset names in metric file (a metadata input file that contains coordinates of nodes).
An example HDF5 file already preprocessed with one HDF5 dataset per trace and day , for the julian day 121 (2018, May, 1st), named ZO_2018_121.h5, could be directly downloaded in the Zenodo deposit associated to our paper.
> curl --head -s https://zenodo.org/record/5645545/files/ZO_2018_121.h5 | grep Content-Length
Content-Length: 7828027009 # file size is ~ 8GB
> wget https://zenodo.org/record/5645545/files/ZO_2018_121.h5
3. Get beamforming_optim code, dependancies installation, compilation and linking
Now that one has the dataset in hdf5 format (it could be our data or any other dataset), we could start the MFP process
3.a. Get beamforming_optim code
The first step is to get the source code. This code is a Fortran code which computes MFP and directly look for local maximum as described in our paper.
> git clone https://gricad-gitlab.univ-grenoble-alpes.fr/lecoinal/resolve.git
The main code for MFP computation with gradient-based optimization is resolve/RESOLVE_beam_src/srcFortran/beamforming_optim.f90.
There are two other main codes in the repository :
- beamforming.f90 : a previous version which computes beamformer on a regular grid (no longer maintained)
- beamforming_simulated_annealing.f90 : another version which computes and optimizes beamformer but with simulated annealing optimization procedure, developed for another project : thèse C. Gradon, Localization of subsurface seismic events on the San Jacinto fault using a dense array of sensors and Gradon et al., 2019).
3.b. Install dependancies
----- option 1
You can install your own dependancies:
For our code to work, one should install the following dependencies:
- gfortran
- hdf5
- lapack
- blas (BLAS/Lapack are only used of case of
ln_csdm=1CSDM built into memory) - nlopt (not used in this paper, we used ASA047 library instead)
One should also install ASA047, which is a FORTRAN77 library which seeks to minimize a scalar function of several variables using the Nelder-Mead algorithm, by R O Neill. ASA047 is Applied Statistics Algorithm 47.
Here is a way to locally install ASA047 library:
# Get src files from https://people.sc.fsu.edu/~jburkardt/f77_src/asa047/asa047.f
# https://people.sc.fsu.edu/~jburkardt/f77_src/asa047/asa047.html
> gfortran -c asa047/asa047.f
> ar qc ${HOME}/softs/asa047/lib/libasa047.a *.o
----- option 2
Or you can use our .yml conda environment and create the MFP conda environment as:
conda env create -f MFP_env.yml
you should still install the asa047 library as:
cd ${HOME}/softs/
wget https://people.sc.fsu.edu/~jburkardt/f77_src/asa047/asa047.f
gfortran -c asa047.f
ar qc libasa047.a *.o
check the path in the Makefile to be sure that you point to the right asa047
3.c. Code compilation and linking
Now, you can compile the code using the Makefile:
> make beamforming_optim
The output should be
h5fc -O2 -cpp -c -o glovarmodule.o ../RESOLVE_beam_src/srcFortran/glovarmodule.f90
h5fc -O2 -cpp -c -o modcomputebeam.o ../RESOLVE_beam_src/srcFortran/modcomputebeam.f90
h5fc -O2 -cpp -c -o modcomputebeam_optim.o ../RESOLVE_beam_src/srcFortran/modcomputebeam_optim.f90
h5fc -O2 -cpp -o beamforming_optim glovarmodule.o modcomputebeam.o modcomputebeam_optim.o ../RESOLVE_beam_src/srcFortran/beamforming_optim.f90 -llapack -lblas -I/usr/include -lnlopt -L/home/lecoinal/softs/asa047/lib -lasa047
3.d. For GriCAD users (Université Grenoble Alpes)
An installation procedure with using Guix environment on GriCAD clusters (only luke and dahu) is detailed here:
Dependencies list, and installation/compilation instructions on GriCAD clusters (guix environment).
4. Prepare metadata
In this step, you should prepare the metadata. In our case, we projected station coordinates on lambertII.
We fix the centre of the array as the 2D-barycentre among the 98 actives stations (indeed the stations AR003 and AR077 never worked and their coordinates were not taken into account to define the center of the array). Our (0,0) relative coordinate corresponds to this UTM coordinate (343058.45425510209, 5091922.07665306050).
We get a coordinates file in relative meters to array center. We keep absolute value for station elevations. You can find our metadata here: Meters Coordinates file for Argentiere ZO_2018 dataset
Here is a map of nodes locations on this relative meter grid.
This can be done with our code
5. Edit the namelist
The namelist is the input file for our code. It contains all necessary information for the code to run (e.g., starting points, 3D space, input data reading management, frequency band, optimization settings, ...).
Here is an example of namelist for computing source location :
- for the first five 1-sec time window with 0.5-sec time-step chunks : i.e. [0:1sec], [0.5:1.5sec], [1:2sec], [1.5:2.5sec], [2:3sec] (see namparam namelist section)
- with coherent (along subfrequencies) bartlett computation (
ln_coh=1) (namcoh namelist section) - with providing (x=20m y=-20m z=2113.73m and c=1800m/s) starting point location (namxy, namHsource and namspeed namelist sections)
- for the central Frequency 17Hz (Fmin - Fmax = 15 - 19 Hz with 0.1Hz subfrequencies : Nf=41) (namindice_freq namelist section)
- with using N=4 parameters (x,y,z,c) and ASA047 library for the optimization (namoptim namelist section)
- with applying a penalty on objective function when optimization path goes above the surface of the glacier (nampenalty namelist section)
Important note: One execution with one namelist corresponds to only one starting point and only one central frequency. The code only iterates over time (here we are doing 5 iterations). If you want to run with different starting points and / or different central frequency, you should run the code several times with different namelist. Each namelist has a set of parameters. We did it this way so the computation can be parallelized with grid computing technique (as embarrassingly parallel jobs).
To give of order of magnitude, our MFP code takes less than 1 hour on 1 core to process 1 hour of data (98 sensors, 500 Hz) with 1-sec time-width and 0.5-sec time-step (i.e. 7200 iterations), for 1 starting point and 1 frequency range with Nf=41 subfrequencies.
Here we have a small number of sensors but the code has already been used on denser networks (e.g. the San Jacinto array with 1108 sensors, or the Glenhaven experiment with 10050 sensors).
When computing MFP in regular grid mode (beamforming.f90) or optimization mode (beamforming_optim.f90 as in our paper), with using ln_coh=1 and ln_csdm=0 options, we measured ~ 140 micro-seconds elapsed to evaluate one MFP value with a configuration with Nr=98 sensors and Nf=41 subfrequencies:
When computing the first 7200 time-steps from the julian day=121:
- in regular grid mode with an (x,y,z,c) regular grid with 3240 grid nodes (Nx=9,Ny=9,Nz=10,Nc=4):
ElapsedTime: 3268.091749242 sec, corresponding to3268 / 7200 / 3240 = 140e-6 sec per MFP evaluation - in optimization mode (with a maximum value of 3000 function evaluations as in our paper):
ElapsedTime: 1986.884179260 sec, and we get a total number of function evaluationssum(infooptim(1:7200,1)) = 14013563, this corresponds to1986 / 14013563 = 141e-6 sec per MFP evaluation
We thus would have a similar CPU cost by running the MFP code with optimization (with maxFunctionEvals=3000) with 29 starting points or by running regular grid MFP code on a regular grid with only ~ 3000x29 grid nodes.
Note about ln_coh parameter:
The code is able to compute incoherent MFP (beamforming incoherent) or coherent MFP (beamforming coherent), depending on how we deal with subfrequencies (see RESOLVE project wiki):
- beamforming incohérent
Bartlett: B = {{1}\over{N_f * N_r^2}} * \Sigma_{f=1}^{N_f} | \Sigma_{i=1}^{N_r} exp^{j(\theta_{data}(i,f)-\psi_{model}(i,f))} |^2
- beamforming cohérent
Bartlett_{coherent}: B_{coh} = {{1}\over{N_f^2 * N_r^2}} * | \Sigma_{f=1}^{N_f} \Sigma_{i=1}^{N_r} exp^{j((\theta_{data}(i,f)-\theta_{data}(iref,f))-(\psi_{model}(i,f)-\psi_{model}(iref,f)))} |^2
Nf is the number of subfrequencies (Nf=41 in our paper)
Nr is the number of nodes (Nr=98 in our paper)
Computing coherent MFP means that CSDM dimensions would be [NrxNf, NrxNf].
In the case of coherent MFP, we choose as iref the node which is the closest to the starting point for the optimization (we choose the minimal horizontal distance, not taking into account node elevations).
In our paper, we always compute coherent Bartlett over subfrequencies.
6. Run the code
Once you have completed the namelist, you can use the wrapper submit.sh to run the MFP code. Using this submission script, you can directly change the starting points and the input/output directory without changing the namelist (you can use namelist as a template that is automatically filled by submission script).
> ./submit.sh
Here is an example of log stdout:
ZO_2018_121.h5 20 -20 2113.73
98 zreseau.txt
freq 15.000000000000000 15.100000000000000 15.199999999999999 15.300000000000001 15.400000000000000 15.500000000000000 15.600000000000000 15.699999999999999 15.800000000000001 15.900000000000000 16.000000000000000 16.100000000000001 16.199999999999999 16.300000000000001 16.399999999999999 16.500000000000000 16.600000000000001 16.699999999999999 16.800000000000001 16.899999999999999 17.000000000000000 17.100000000000001 17.199999999999999 17.300000000000001 17.399999999999999 17.500000000000000 17.600000000000001 17.699999999999999 17.800000000000001 17.899999999999999 18.000000000000000 18.100000000000001 18.199999999999999 18.300000000000001 18.399999999999999 18.500000000000000 18.600000000000001 18.699999999999999 18.800000000000001 18.899999999999999 19.000000000000000
Read the station list and their geographical location from zreseau.txt
Read input data: offset,count: 0 1500
1-based index of reference sensor is 44
and its horizontal distance from starting point is 24.711173707454687 m
1-based index ref_inode_coh 44 ZO.AR045.00.DPZ.D.2018.121
read inputs and init: 9.76562500E-02
it | totsec : 5 1.76171875
write outputs: 0.00000000
Write_output 20211111 103231.800 0.00000000
ElapsedTime: 1.891635726 sec
Now the MFP process is finished. You can re-run it with different parameters.
For instance, we try to compute MFP for one hour of data between 14:00:00 and 15:00:00 on April, 25th, 2018 (julday=115).
We thus change namelist as :
[...]
jbeg = 25200000 ! 0 ! The starting point (in nb of samples, not in sec) (Caution: index starts at 0)
[...]
nit = 7200 ! 5 ! 43200 ! 5 ! 172799 ! 43200 ! 18000 ! 86399
/
[...]
/
&namio
nn_io = 7200 !
1-hour with 0.5 sec timestep = 7200 time-steps
jbeg = 14 x 3600 x 500Hz = 25,200,000
and change input file name in submit.sh:
IFILE="ZO_2018_115.h5"
And submit again
> ./submit.sh
ZO_2018_115.h5 20 -20 2113.73
98 zreseau.txt
Read the station list and their geographical location from zreseau.txt
Read input data: offset,count: 25200000 1800250
1-based index of reference sensor is 44
and its horizontal distance from starting point is 24.711173707454687 m
1-based index ref_inode_coh 44 ZO.AR045.00.DPZ.D.2018.115
read inputs and init: 3.63671875
it | totsec : 500 137.046875
it | totsec : 1000 137.910156
it | totsec : 1500 135.195312
it | totsec : 2000 133.179688
it | totsec : 2500 141.476562
it | totsec : 3000 134.015625
it | totsec : 3500 131.273438
it | totsec : 4000 140.906250
it | totsec : 4500 140.867188
it | totsec : 5000 159.601562
it | totsec : 5500 166.832031
it | totsec : 6000 153.167969
it | totsec : 6500 159.898438
it | totsec : 7000 155.480469
it | totsec : 7200 65.3828125
write outputs: 3.90625000E-03
Write_output 20211103 160718.892 3.90625000E-03
No all objects to modify layout
All objects to apply filter are...
All with GZIP, parameter 6
Making new file ...
-----------------------------------------
Type Filter (Compression) Name
-----------------------------------------
group /
dset GZIP (1.115:1) /Baval
dset GZIP (3.897:1) /infooptim
dset GZIP (1022.238:1) /patchesInd
dset GZIP (1.087:1) /xyzc
ElapsedTime: 2095.939036131 sec
7. Investigate the outputs
The output of the MFP process are in .h5. You can explore the output using:
> h5ls -d out_ZO_2018_121.h5
We gives
- The values of the MFP output (the best phase coherence), here you have 5 localizations corresponding to the 5 iterations in time.
Baval Dataset {5}
Data:
(0) 0.00705546, 0.0103456, 0.00935663, 0.0118354, 0.003344
- The information about the optimization. First column is the number of function evaluations before the optimization algorithm converges. If the iteration number is larger than 3000 (this stopping condition could be changed in the namelist), this means that the optimization did not converged. Second column is the number of algorithm re-start. Third column is the status code of the optimization algorithm library (see ASA047 documentation), for example a value of 2 for the status code means that we reach the stopping condition for maximum number of function evaluations (example with first, second and last iterations here). Third column of infooptim could thus be used as a mask (if status code different from zero) on the MFP outputs.
infooptim Dataset {3, 5}
Data:
(0,0) 3021, 3041, 1755, 1776, 3074, 17, 17, 9, 9, 17, 2, 2, 0, 0, 2
- The information about the availability of input data. Remember that we preprocessed input data by filling gaps with zero. Here we are computing MFP with 1-sec time width. If, for the current iteration, one node has only zero value during the whole 1-sec time-width iteration, we store this information here. Thus patchesInd would be a binary mask that we could apply on the MFP outputs.
patchesInd Dataset {5, 98}
Data:
(0,0) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
(0,40) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
(0,80) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
(1,22) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
(1,62) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
(2,4) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
(2,44) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
(2,84) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
(3,26) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
(3,66) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
(4,8) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
(4,48) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
(4,88) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
- The localization (x, y, z) of the sources, as well as the associated phase velocity.
xyzc Dataset {4, 5}
Data:
(0,0) 27.8916, -36.5648, -29.923, -7.94212, 53.7745, 31.1171, 32.0942, 89.3524, 6.01179, -28.3819, 978.376, 1465.83, 1972.66,
(2,3) 2366.07, 1597.28, 463.052, 932.897, 2596.6, 1554.65, 458.929
The complete output of our investigation (29 starting point, 16 Frequency, 34 days, 98 sensors) can be found on our Zenodo deposit, together with the associated Matlab codes to read the .h5 files, plot the data and create source location density maps as in our paper.
8. Big data computation using grid computing technique
The code is purely sequential. One execution of the code iterates over the time-steps inside one single input file (daily input files in this example). Running several frequency bands and/or several starting points for the optimization and/or several input files (several days of input data in this example) means executing several times the same code with different namelist parameters.
To do that (sensitivity analysis through parametric study along frequency bands and starting points, combined with massive data processing along days of experiment), we used a grid computing technique. The middleware CiGri allows us to launch bags of embarrassingly jobs to compute MFP for 34 days, 16 frequency bands, 29 starting points = 15776 grid jobs.
We provide a CiGri parameter file: it is text file with one line of parameters corresponding to one CiGri job associated to these parameters. This parameters line becomes the arguments ($1 $2 $3 ...) of the bash submission script. CiGri submits jobs on targeted clusters. The bash submission script updates the namelist.skel template with replacing patterns in the namelist by the parameters line. Here we use this syntax <<<IFILE>>> or <<<XSTART>>> ... in namelist.skel template to indicate that this pattern will be updated by submission script.
The parameters (namelist and CiGri parameters file) and the source code revision used in our paper are detailed in the wiki corresponding to the computations performed for this paper CiGri experiment number 15423.
The nn_io parameter from namelist is useful in the case of a massive deployment of the code over a grid computing infrastructure. This parameter allows to set the pre-load of input data in case of several time-widths iterations. The main part of memory footprint of code execution is associated to the input read (here 98 nodes with 500Hz sampling rate), in particular as we use the code with setting ln_csdm=0 (we don't build CSDM to compute Bartlett, we use the simple formulation with
bartlett = z*conj(z) with z = sum (exp (j*(phi_data-phi_model(:))))
as we don't care about the amplitudes, we only consider phases and thus we only allocate 2 vectors with nbnodes = 98 elements or 98 x 41 subfrequencies of data phases phi_data and replica phases phi_model). Reading input file is thus conditioning the behavior of a job in terms of IOps and memory footprint. Setting the parameter nn_io=1 means that we read input file at each iteration (here one iteration is a 1-sec time-width with a 0.5-sec time-step). nn_io=1 is therefore not optimal in this context where we overlap the input data by 50%. Besides, in case of massive grid computing deployment (if a large amount of grid jobs are running simultaneously because of one or several thousands computing cores are free), a too low nn_io value can lead to shared filesystem IOps saturation and slow down the file system for all users. On the contrary, a too high nn_io value would lead to a memory overload for the grid job. Indeed, reading a full input daily file at first iteration would lead to a memory footprint around 16GB (one input daily file is 8GB but HDF5 datasets are 200% gzip compressed : 98nodes * 86400sec * 500Hz * 4 (float32) = 15.8GB). In the grid clusters from GriCAD infrastructure, that we used in our paper, the computing nodes are heterogeneous among the 4 grid clusters (that is a feature inherent in grid computing, be able to take advantage of a heterogeneous clusters infrastructure), the cores with the least amount of by-core memory have 2GB memory per core. We thus set the nn_io value to find a compromise to limit the IOps on the shared filesystem while limiting the memory footprint of an execution below 2GB. In our computations for the paper, we used nn_io=7200 (which corresponds to 1-h input data pre-load, with 0.5-sec time-step), which resulted in a grid job consuming 700MB RAM (98nodes * 500Hz * 3600sec * 4(F32) = 700MB).
For these computations, whose runtime depends on convergence speed for the optimization, grid computing technique is optimal in the sense of synchronization between jobs is not required. Each job is totally independant from others. A more classical parallelization technique (such as MPI) would be less effective to implement because of the imbalance between MPI processes inherent to our jobs.
For our paper, the CiGri experiment contains 15776 grid jobs. Each job was a sequential monocore job, with an average duration (elapsed time) of 12 hours per job (depending on the convergence speed for the optimization). We targeted 4 different clusters available through CiGri, with heterogeneous resources. We compute best-effort type jobs : that means that our grid jobs could be killed at any time by other jobs with more priority (because of their owner's privileges related to the financing of clusters opened to cigri, or because the other jobs are massively parallel). Our grid campaign was 7 days elapsed and used a total of 190,000 CPU-hours (core-hours).
To recombine the 15776 outputs from the 15776 grid jobs, we used Python script and h5py module. Recombining 15776 outputs to create one single ~50GB output file is a sequential monocore ~7h elapsed job.