Difference between revisions of "Computing a theoretical dispersion curve"

From GeopsyWiki
Jump to navigation Jump to search
 
(15 intermediate revisions by 2 users not shown)
Line 1: Line 1:
 
[[Category:Tutorials]]
 
[[Category:Tutorials]]
  
This tutorial details all the steps needed to compute a dispersion with '''gpdc'''.
+
This tutorial details all the steps needed to compute a dispersion with '''gpdc'''. Another more interactive way to compute dispersion curves is [[Gplivemodel]].
  
 
== Creating a model file ==
 
== Creating a model file ==
  
Crete a text file (e.g. "test.model"):
+
Crete a text file (e.g. "[[Media:Test.model|test.model]]"):
  
 
  # My first model: two layers over a half-space
 
  # My first model: two layers over a half-space
Line 13: Line 13:
 
  # Thickness(m), Vp (m/s), Vs (m/s) and density (kg/m3)
 
  # Thickness(m), Vp (m/s), Vs (m/s) and density (kg/m3)
 
  7.5  500  200 1700
 
  7.5  500  200 1700
  25  1350  190 1900
+
  25  1350  210 1900
 
  # Last line is the half-space, its thickness is ignored but the first column is still mandatory
 
  # Last line is the half-space, its thickness is ignored but the first column is still mandatory
 
  0  2000 1000 2500
 
  0  2000 1000 2500
Line 23: Line 23:
 
To get the dispersion curve of Rayleigh fundamental mode:
 
To get the dispersion curve of Rayleigh fundamental mode:
  
  gpdc test.model
+
  gpdc Test.model
  
 
You get as output:
 
You get as output:
Line 35: Line 35:
 
  # CPU Time = 1 ms
 
  # CPU Time = 1 ms
 
  # Mode 0
 
  # Mode 0
  0.2 0.00107882203003697
+
  0.2 0.00107875907546066
  0.209523150557933 0.00107913814934707
+
  0.209523150557933 0.00107907141898638
 
  [...]
 
  [...]
  19.0909691332367 0.00533451456923106
+
  19.0909691332367 0.00525624231593241
  20 0.00533094289276869
+
  20 0.00526288933387501
  
 
To compute Rayleigh higher modes (4 in this case) as well as the fundamental mode:
 
To compute Rayleigh higher modes (4 in this case) as well as the fundamental mode:
  
  gpdc test.model -R 5
+
  gpdc Test.model -R 5
  
 
To compute Love higher modes (4 in this case) as well as the fundamental mode:
 
To compute Love higher modes (4 in this case) as well as the fundamental mode:
  
  gpdc test.model -R 0 -L 5
+
  gpdc Test.model -R 0 -L 5
  
To compute Love modes, you must set explicitly the number of Rayleigh to 0, it is one by default.  
+
To compute Love modes, you must set explicitly the number of Rayleigh to 0, it is one by default.
  
 
== Saving the dispersion curves ==
 
== Saving the dispersion curves ==
  
  gpdc test.model -R 2 > test_Rayeigh_2modes.disp
+
  gpdc Test.model -R 2 > Test_Rayeigh_2modes.disp
  
 
== Plotting the dispersion curves ==
 
== Plotting the dispersion curves ==
  
We use '''figue''' to plot couples of values (X, Y), additionally we provide a '''make-up''' file to get the correct type and labeling of axes. Before running the next command you can download [dc.mkup].
+
[[Image:gpdc_dispersion_plot.png|thumb|right|300px|Dispersion curve plot with figue]]
  
gpdc test.model -R 2 | figue -c -m dc.mkup
+
We use '''figue''' to plot couples of values (X, Y), additionally we provide a '''make-up''' file to get the correct type and labeling of axes. Before running the next command you can download [[Media:Dc.mkup|Dc.mkup]].
  
[[Image:gpdc_dispersion_plot.png|thumb|right|400px|Dispersion curve plot with figue]]
+
gpdc Test.model -R 2 | figue -c -m Dc.mkup
  
 
== Going further ==
 
== Going further ==
Line 73: Line 73:
 
=== Programming interface ===
 
=== Programming interface ===
  
The computation of dispersion is implemented in ''libqtbwave.so'' (or equivalent for Windows and Mac). You can access it through various languages: C++, C or Fortran
+
The computation of dispersion is implemented in ''libQGpCoreWave.so'' (or equivalent for Windows and Mac). You can access it through various languages: C++, C or Fortran. Complete examples are provided for C and Fortran in the source distribution archive. Download ''gpdc'' to have the simplest archive.
  
A C++ example:
+
==== A C++ example ====
  
  #include "qtblayeredmodel.h"
+
  #include <QGpCoreWave.h>
#include "qtbrayleigh.h"
 
//#include "qtblove.h"
 
#include "qtbdispersion.h"
 
 
  [...]
 
  [...]
 
  // 3-layer model initialization
 
  // 3-layer model initialization
  QtbLayeredModel model(3);
+
  LayeredModel model(3);
 
  model.h()[0]=7.5;
 
  model.h()[0]=7.5;
 
  model.h()[1]=25.0;
 
  model.h()[1]=25.0;
Line 101: Line 98:
 
  [...]
 
  [...]
 
  x[49]=...;
 
  x[49]=...;
  // Or use QtbCurve to generate linear or log series
+
  // Or use Curve to generate linear or log series
  QtbCurve<QtbPoint1D> c;
+
  Curve<Point1D> c;
 
  c.line( minFreq, 0.0, maxFreq, 0.0 );
 
  c.line( minFreq, 0.0, maxFreq, 0.0 );
  c.resample( nSamples, minFreq, maxFreq, Qtb::Log | Qtb::Function );
+
  c.resample( nSamples, minFreq, maxFreq, Log | Function );
 
  c.multiply( 2*M_PI ); // convert to angular frequency
 
  c.multiply( 2*M_PI ); // convert to angular frequency
 
  QVector<double> x = c.xVector();
 
  QVector<double> x = c.xVector();
 
  // Dispersion computation, fundamental mode only
 
  // Dispersion computation, fundamental mode only
  QtbRayleigh rayleigh( &model );
+
  Rayleigh rayleigh( &model );
  //QtbLove love( &model );
+
  //Love love( &model );
  QtbDispersion dispersion( 1, &x);
+
  Dispersion dispersion( 1, &x);
 
  dispersion.calculate( modelRayleigh, 0 ); // 0 is to avoid ellipticity computation
 
  dispersion.calculate( modelRayleigh, 0 ); // 0 is to avoid ellipticity computation
 
  // Access computed values
 
  // Access computed values
  QtbCurve<QtbPoint2D> * dc = dispersion.curve(0);
+
  Curve<Point2D> * dc = dispersion.curve(0);
  // See QtbCurve header for more information
+
  // See Curve header for more information
 
  // Or directly, fundamental mode, ith sample value.
 
  // Or directly, fundamental mode, ith sample value.
 
  dispersion.mode(0)[i]
 
  dispersion.mode(0)[i]
Line 120: Line 117:
 
To compile:
 
To compile:
  
  g++ [...] -lqtbtools -lqtbwave
+
  g++ [...] -lQGpCoreTools -lQGpCoreWave
  
A Fortran example:
+
==== Fortran interface ====
  
  call dispersion_curve_init
+
  call dispersion_curve_init(verbose)
  call dispersion_curve_rayleigh(nLayers, h, vp, vs, rho, nSamples, omega, nModes, slowness );
+
  call dispersion_curve_rayleigh(nLayers, h, vp, vs, rho, nSamples, omega, nModes, slowness, group);
  call dispersion_curve_love(nLayers, h, vp, vs, rho, nSamples, omega, nModes, slowness );
+
  call dispersion_curve_love(nLayers, h, vp, vs, rho, nSamples, omega, nModes, slowness, group);
  
 
{| border="1" cellpadding="5" cellspacing="0"
 
{| border="1" cellpadding="5" cellspacing="0"
Line 133: Line 130:
 
  ! Type
 
  ! Type
 
  ! Description
 
  ! Description
 +
|-
 +
| verbose
 +
| Integer 4 bytes
 +
| 0 minimal ouput, 1 verbose output
 
  |-
 
  |-
 
  | nLayers
 
  | nLayers
Line 169: Line 170:
 
  | Vector of floats 8 bytes (nSamples*nModes, elements)
 
  | Vector of floats 8 bytes (nSamples*nModes, elements)
 
  | Output of slowness values
 
  | Output of slowness values
 +
|-
 +
| group
 +
| Integer
 +
| 0 for phase, 1 for group
 
  |}
 
  |}
  
Line 174: Line 179:
 
To compile:
 
To compile:
  
  gfortran [...] -lqtbtools -lqtbwave
+
  gfortran [...] -lQGpCoreTools -lQGpCoreWave
  
 
=== Theoretical considerations ===
 
=== Theoretical considerations ===
  
The equations and methods implemented in ''libqtbwave.so'' are fully described from a theoretical point of view in [http://marc.geopsy.org/publications.html Marc Wathelet's PhD thesis] (chapter 3).
+
The equations and methods implemented in ''libQGpCoreWave.so'' are fully described from a theoretical point of view in [http://marc.geopsy.org/publications.html Marc Wathelet's PhD thesis] (chapter 3).

Latest revision as of 14:41, 11 January 2013


This tutorial details all the steps needed to compute a dispersion with gpdc. Another more interactive way to compute dispersion curves is Gplivemodel.

Creating a model file

Crete a text file (e.g. "test.model"):

# My first model: two layers over a half-space
# First line: number of layers
3
# One line per layer:
# Thickness(m), Vp (m/s), Vs (m/s) and density (kg/m3)
7.5  500  200 1700
25  1350  210 1900
# Last line is the half-space, its thickness is ignored but the first column is still mandatory
0   2000 1000 2500

All lines beginning with '#' are considered as comments and they are ignored. Several models can be chained in one single file, but no comment line is allowed between layers. Column separators can be either TAB or SPACE. The number of consecutive separators does not matter.

Computing the dispersion curves

To get the dispersion curve of Rayleigh fundamental mode:

gpdc Test.model

You get as output:

# My first model: two layers over a half-space
# First line: number of layers
# One line per layer:
# Thickness(m), Vp (m/s), Vs (m/s) and density (kg/m3)
# Last line is the half-space, its thickness is ignored but the first column is still mandatory
# 1 Rayleigh dispersion mode(s)
# CPU Time = 1 ms
# Mode 0
0.2 0.00107875907546066
0.209523150557933 0.00107907141898638
[...]
19.0909691332367 0.00525624231593241
20 0.00526288933387501

To compute Rayleigh higher modes (4 in this case) as well as the fundamental mode:

gpdc Test.model -R 5

To compute Love higher modes (4 in this case) as well as the fundamental mode:

gpdc Test.model -R 0 -L 5

To compute Love modes, you must set explicitly the number of Rayleigh to 0, it is one by default.

Saving the dispersion curves

gpdc Test.model -R 2 > Test_Rayeigh_2modes.disp

Plotting the dispersion curves

Dispersion curve plot with figue

We use figue to plot couples of values (X, Y), additionally we provide a make-up file to get the correct type and labeling of axes. Before running the next command you can download Dc.mkup.

gpdc Test.model -R 2 | figue -c -m Dc.mkup

Going further

Help about command line options

To get help about all options:

gpdc -h all

Programming interface

The computation of dispersion is implemented in libQGpCoreWave.so (or equivalent for Windows and Mac). You can access it through various languages: C++, C or Fortran. Complete examples are provided for C and Fortran in the source distribution archive. Download gpdc to have the simplest archive.

A C++ example

#include <QGpCoreWave.h>
[...]
// 3-layer model initialization
LayeredModel model(3);
model.h()[0]=7.5;
model.h()[1]=25.0;
model.slowP()[0]=1.0/500.0;
model.slowP()[1]=1.0/1350.0;
model.slowP()[2]=1.0/2000.0;
model.slowS()[0]=1.0/200.0;
model.slowS()[1]=1.0/210.0;
model.slowS()[2]=1.0/1000.0;
model.rho()[0]=1800;
model.rho()[1]=1900;
model.rho()[2]=2500;
// Initialize frequency sampling, in angular frequency (omega)!!
QVector<double> x;
x.resize(50);
x[0]=...;
[...]
x[49]=...;
// Or use Curve to generate linear or log series
Curve<Point1D> c;
c.line( minFreq, 0.0, maxFreq, 0.0 );
c.resample( nSamples, minFreq, maxFreq, Log | Function );
c.multiply( 2*M_PI ); // convert to angular frequency
QVector<double> x = c.xVector();
// Dispersion computation, fundamental mode only
Rayleigh rayleigh( &model );
//Love love( &model );
Dispersion dispersion( 1, &x);
dispersion.calculate( modelRayleigh, 0 ); // 0 is to avoid ellipticity computation
// Access computed values
Curve<Point2D> * dc = dispersion.curve(0);
// See Curve header for more information
// Or directly, fundamental mode, ith sample value.
dispersion.mode(0)[i]

To compile:

g++ [...] -lQGpCoreTools -lQGpCoreWave

Fortran interface

call dispersion_curve_init(verbose)
call dispersion_curve_rayleigh(nLayers, h, vp, vs, rho, nSamples, omega, nModes, slowness, group);
call dispersion_curve_love(nLayers, h, vp, vs, rho, nSamples, omega, nModes, slowness, group);
dispersion_curve_rayleigh arguments
Argument Type Description
verbose Integer 4 bytes 0 minimal ouput, 1 verbose output
nLayers Integer 4 bytes Number of layers
h Vector of floats 8 bytes (nLayers elements) Thicknesses of layers (m)
vp Vector of floats 8 bytes (nLayers elements) Vp in each layer (m/s)
vs Vector of floats 8 bytes (nLayers elements) Vs in each layer (m/s)
rho Vector of floats 8 bytes (nLayers elements) Density in each layer (kg/m3)
nSamples Integer 4 bytes Number of frequency samples
omega Vector of floats 8 bytes (nSamples elements) Angular frequencies (rad/s)
nModes, Integer 4 bytes Number of modes including fundamental
slowness Vector of floats 8 bytes (nSamples*nModes, elements) Output of slowness values
group Integer 0 for phase, 1 for group


To compile:

gfortran [...] -lQGpCoreTools -lQGpCoreWave

Theoretical considerations

The equations and methods implemented in libQGpCoreWave.so are fully described from a theoretical point of view in Marc Wathelet's PhD thesis (chapter 3).