This web pages describes how to use
the software related to the paper by Debayle & Sambridge (2004).
These software are available on request from the first author Eric Debayle. A tar
file will be provided with DS2004 as root
directory. We are interested in obtaining feedback on people's
experiences with our code especially on the performance and behavior of our tomography code regiostack6_ompV3d.f. Any suggestion or improvements are welcome. However, the purpose
of
this web page is to reduce the likelihood of running into problems. We
would be grateful if you do not contact the first author with queries
on the installation, or running of the code unless you have thoroughly
looked for the answer yourself (note that either these HTML pages, the
example files and the source code itself may provide an answer).
This page will be updated from time to time. Last
update is 28th January 2005
Reference
We are happy to make our codes available,
hoping that people give proper reference to our work :
-
Debayle E. and M. Sambridge, Inversion of
massive surface wave data sets : Model construction and Resolution
assessment, Journal of Geophysical Research,
109, B02316,
doi:10.1029/2003JB002652 , 2004. (pdf)
Tomography code for
massive surface wave data sets
Using our tomography code is very simple.
It allows to invert path average phase or group velocities measurements
at different periods or path average shear velocity models at different
depths. The inversion can be performed for the local velocity alone or
for the local velocity plus the 2theta anisotropic terms that are
relevant for Rayleigh and SV waves.
The code can be used indifferently for regional or global tomography.
You don't have to enter the boundaries of the regions you are working
on : because our code explores an "influence zone" around each ray
path, computations are automatically restricted
in the region where the rays are. An output model will always be
obtained for the whole Earth. In parts of the Earth that are not
sampled by any rays the output model stays on the a priori model (i.e.
no lateral variations). For a "regional tomography" you get a "whole
Earth" model with no lateral variations in regions not sampled. You may
therefore just plot the region you are interested in.
The same code can be used either on a single processor PC or on a
multi-processor machine under OpenMP. In this case you just need to
compile it with the OpenMP options (see example for an IBM Power 4
computer in the compilation section). The
performances should be optimal with 2 threads.
SOURCE CODES
The source codes for tomography are in DS2004/TOMO/src.
INPUT FILES
2 input files are needed by regiostack6_ompV3d.f :
A typical
'intomodesVs' file with 3 paths has the following format :
-------------------------------------------------------
z1 path 1 "z"=vertical
"1": number of path; "path 1"=path name
-8.4600 107.9600 21.4233 -158.0150 Lat lon event, lat lon station
4.4236 4.4115 4.3850... velocities
for each layers 1, 2, 3...
0.0389 0.0389 0.0389... error
on velocities for each layers 1, 2, 3...
z2 path 2
-8.2927 108.4918 -20.0882 146.2545
4.3782 4.4208 4.4725...
0.0237 0.0268 0.0260...
z3 path 3
-10.5104 118.7370 -32.9266 117.2333
4.5507 4.5817 4.6191...
0.0298 0.0223 0.0207...
------------------------------------------------------
Structure of an "I.vs"
file with 3 paths :
---------------------------------------------------
0.00 Additional error that you may want
to add on your data (*)
2.00 Step in degree along paths
3 Number of paths (correspond to
parameter NTRA in the code)
2.00 Size of the cell for output vsfin and
anvsfin files in degree
3 Nber total of depths/periods in the
intomodesVs file
2 1 2 Nber and indice of depths/periods to
invert (here only the velocities at 40. and 50. km are inverted)
40. 50. 75. Depths/periods of all the
layers present in the intomodesVs file
intomodesVs Name of input intomodesVs file
vsfin Output file for velocity
anvsfin Output file for 2theta
anisotropic parameter A1,A2
400. 400. Correlation length (km) for
isotropic and 2theta anisotropic parameters
40.00 4.327 .05 0. 0. 0.005
Depths/periods,
model a priori velocity, model a priori standard deviation for
velocity, a priori values for model parameters A1, A2 (0. 0. means no a
priori anisotropy), model a priori standard deviation for the 2theta
anisotropic parameters A1 and A2 (0. would mean that you perform an
isotropic inversion).
50.00 4.323 .05 0. 0. 0.005
75.00 4.329 .05 0. 0. 0.005
---------------------------------------------------
(*) You may think
that some uncertainties in the origin time, in the source mechanism or possible
mislocations in the epicenter may not be reflected in the path average
errors provided in the
intomodesVs file (this may be the case with the a posteriori error
issued from a waveform inversion).
In this case you have the option to add an error in km/s that will be
added to the data error in
each path assuming the errors are normally distributed.
COMPILATION
The source codes for tomography are in DS2004/TOMO/src.
-
" make regiostack6_ompV3d.ifc"
to compile with the ifc compiler (ttp7 is a flag that optimizes for Pentium 4 processor and may need to be edited to
suite the users platforms).
-
" make regiostack6_ompV3d.f90 " to compile with the
pgf90 compiler.
For multiprocessors applications you need
to link with the OpenMP libraries of your
compiler. An example is given in Makefile_omp for the IBM compiler used
on the IBM power 4 supercomputer used at IDRIS (Institut du Développement et des ressources
en informatique scientifique). Once compiled with OMP the performances
of the code should be optimal with 2 threads.
RUN AN EXAMPLE
An example of input files is provided in directory DS2004/TOMO/EXAMPLE.tomo. Try
:
The inversion of this synthetic dataset of 5672 paths
average S velocity models calculated in the 3SMAC model takes
about 7 min on my Pentium
4
PC 2.4 Ghz processor with 2Gbytes of RAM. For a given number of rays,
size of the blocks, step along paths and correlation length the
computation time of our tomography code depends on the geometry of the
rays coverage. The performances of the code
announced in
our JGR paper are a lower bound and you may find that using an
equivalent computer, in most case our code is faster.
The so file contains various information regarding the results. Among
the interesting ones :
In most case the code should work without any changes.
REPORTED PROBLEMS
3 parameter corresponding to the dimension
of some arrays could not be allocated dynamically :
"NTR must be
increased" : NTR is the maximal number of possible rays in a cell. It
is arbitrary fixed according to NTR which is the total number of rays
in the loop from line 245 to 255 of the code. If you've got a lot of
memory you may just set NTR=NTRA and forget it. If you need to save
memory a technique is to make a guess (as I do in the loop from line
245 to 255 ) or if possible to do a first run and a ''grep Number so
''which provide you with the ''Number max of paths in a cell'' and the
"Number used for arrays allocation" (i.e. the NTR you entered). NTR
will vary depending on ray geometry and block size
"Erreur SA, IJA" : NF is used in arrays SA IJA and allow to store in
compact form the data by data matrix that is inverted. The technique
used is a "row indexed compact sparse matrix storage technique "
described in numerical recipes Fortran 77 chapter 2.7 page 71.
Depending on how sparse is the matrix we can store the data by data
matrix (i.e. NTRA*NTRA) in a (NTRA*NF) vector with NF<< NTRA. I
propose values of NF that work in most case in the loops from line 238
to 242 of the code. If the code stops with "Erreur a L'allocation de
SA" or "Erreur a L'allocation de IJA" it probably means that you don't
have enough memory to allocate SA and IJA so that you should decrease
NF (and eventually NTR after a ''grep Number so '' as described
above)or buy more memory. If you get the message "Size arrays SA IJA
too small value of NF must be increased in main code" try to increase
NF. After a first successful run you can
know
exactly what are the actual size and the allocated size for SA and IJA.
Do a "grep Size so".
"Stack Full error" : Stack_size is used in the LIFO (Last In First Out)
routines to store the cells contained in the "influence zone" of a
given point of the path. The problem occur when you use a very large
horizontal correlation length or a small cell size (or both). In this
case just increase the value of the "stack_size" parameter which is set
to 1000 in the main code, the CmGtstack2, the GcmGtstack, the pushf and
the popf routines.
PLOTTING THE
RESULTS
After running the
example, the output velocities and anisotropy values are stored in the
vsfin and anvsfin files located in DS2004/TOMO/EXAMPLE.tomo.
- 1. Quick plot :
-
- Use regio62gmt.f
and answer the questions (an example is given below)
- --------------------------------------------------------
../bin/regio62gmt
On prends 1066a comme modele de reference?(1=oui)
0
entrer le pas en degres du fichier initial le
nombre de lignes et de colonnes du fichier
2 90 180
entrer le fichier vitesse
vsfin
entrer le fichier aniso
anvsfin
vsfin anvsfin
entrez le nombre de prof des fichier
1
Entrez corr, facteur de correction pour le plot de
l'anisotropie(corr=5->1cm = 5% d'aniso sur les cartes)
5
entrez 1 pour faire des coupes en profondeur
(il faut des fichiers Interpolles)
(ne marche que dans le cas de Vs)
0
Voulez-vous lire la densite et Vs aux prof de 1066a (1)
ou dans PREM smoothe et interpolle tout les 25 km (2)?
2
-------------------------------------------------------
-
- Try "csh
figvsa.gmt 01" to plot the result of the first layer. Click here to see the result
that you should obtain.
N.B. : Note that the result shows how the
isotropic model a priori 3SMAC by Nataf and Ricard (1996) is retrieved
with an arbitrary synthetic coverage of 5672 paths, when anisotropy is
allowed in the inversion. The anisotropy displayed is an artifact
related to a weak azimuthal distribution of rays but this result is
completely dependent of the arbitrary ray
coverage and a priori chosen. The ability of the technique to retrieve
an input model depends on the data sampling and a priori chosen, and it
will be almost perfect if the data coverage is excellent and does not
require strong a priori damping. This example is just a quick way to
try a run of the codes and should not be interpreted.
-
Informations to plot the results
with your own codes
tt = period or depths
pxyo = size of the cell in line 4 of
I.vs
nxo = number of lines
(int(180./pxyo))
nyo = number of columns
(int(360./pxyo))
ulf = Co(T) (equation 18) or
beta_v0(z) (equation 19 of Debayle & Sambridge 2004)
dan1 = A1 (equation 18) or
G_c/(2*rho*beta_v0(z)) (equation 19 of Debayle & Sambridge 2004)
dan2 = A2 (equation 18) or
G_s/(2*rho*beta_v0(z)) (equation 19 of Debayle & Sambridge 2004)
The file vsfin (logical unit 7) contain
the local velocity
and is written according to the following the format :
Loop on period/depths
write(7,*)tt
do 127 i=1,nxo
127 write(7,2014)(ulf(i,j),j=1,nyo)
End loop on period/depths
and anvsfin (logical unit 9) contain the
anisotropic
parameters and is written according to the following the format :
Loop on period/depths
write(9,*)tt
do 125 i=1,nxo
125 write(9,2014)(dan1(i,j),j=1,nyo)
do 124 i=1,nxo
124 write(9,2014)(dan2(i,j),j=1,nyo)
End loop on period/depths
with
2014 format(20(f10.6,1x))
Values are stored from the upper left (North pole) to the lower right
(South pole) and from -180 degres to +180 degres.
2D Voronoi diagrams for
surface wave tomography
The building of an optimized Voronoi
Diagram according to Debayle and Sambridge (2004) is an automated
process that requires 2 input files :
Examples of the input files are provided in the directory DS2004/VORONOI/EXAMPLE.voro. The LISTNODES file is
designed to start from a 2*2 degree grid. The
intomodesVs file contains the ray coverage. The intomodesVs file provided in EXAMPLE.voro
(named
intomo3SMAC) is the same file as the one used in DS2004/TOMO/EXAMPLE.tomo.
CODE DESCRIPTION
The main source codes are located in DS2004/VORONOI/src :
COMPILATION
Go to DS2004/VORONOI/qhull_src. If
you've got the icc compiler try :
-
"icc *.c -c" to build the
object files of the qhull routines geom.c, globals.c, io.c, mem.c,
poly.c, qhull.c, set.c.
N.B. : depending on the compiler
you
may encounter difficulties to compile qhull.c. It has been
working
successfully on sun using
gcc 3.2.1 gcc 2.95.3 and on linux using gcc 3.2.2 (20030222).
In DS2004/VORONOI/src
open the Makefile and set the paths to your Fortran and C compiler.
For a separate
compilation :
RUN AN EXAMPLE AND PLOT THE
RESULTS
After compilation, try
the example in DS2004/VORONOI/EXAMPLE.voro .
"csh
autodelrand3.run intomo3SMAC"
This example build the optimized Voronoi diagram for a synthetic
distribution of 5672 paths. The cshell script autodelrand3.run
produces
two '*.xy' GMT files that contain the starting
(voronoi.init.intomo3SMAC.xy) and optimized (voronoi.intomo3SMAC.xy)
Voronoi diagrams to be plot by GMT. The last line in the cshell script autodelrand3.run call a gmt script vorfin.gmt to plot the results on the screen. vorfin.gmt work if the GMT
software is installed on your computer. It uses ggv to plot the results but this
may be edited to suits the user's platform. Click here to see the results of a successful run
of autodelrand3.run.
top of page