Module rgrid

Raytracing on rectilinear grids

This module contains two classes to perform traveltime computation and raytracing on rectilinear grids:

Grid2d for 2D media
Grid3d for 3D media

Three algorithms are implemented:

the Shortest-Path Method
the Fast-Sweeping Method
the Dynamic Shortest-Path Method

Slowness model can be defined in two ways:

slowness constant within the voxels of the grid (the default)
slowness defined at nodes of the grid

This code is part of ttcr ( https://github.com/groupeLIAMG/ttcr )

ttcrpy.rgrid.Grid2d(x, z, ..., dtype=np.float64) → Grid2d_d or Grid2d_f

Factory that returns a 2D rectilinear-grid raytracer with the requested numeric precision.

Parameters:

x (array-like) – Node coordinates. Converted to the requested dtype automatically.
z (array-like) – Node coordinates. Converted to the requested dtype automatically.
dtype (numpy dtype) – np.float64 (default) for double precision (Grid2d_d), np.float32 for single precision (Grid2d_f).
fsm_gpu (bool) – use OpenCL implementation to run the Fast Sweeping Method on GPU (FSM only, default is False)
class. (All other parameters are forwarded unchanged to the underlying)

Return type:

Grid2d_d or Grid2d_f

class ttcrpy.rgrid.Grid2d_d

class to perform raytracing with 2D rectilinear grids (double precision)

Variables:

x (np.ndarray) – node coordinates along x
z (np.ndarray) – node coordinates along z
dx (float) – node separation along x
dz (float) – node separation along z
shape ((int, int)) – number of parameters along each dimension
nparams (int) – total number of parameters for grid
n_threads (int) – number of threads for raytracing
Constructor
Grid2d (Grid2d(x, z, n_threads=1, cell_slowness=1, method='SPM', aniso='iso', eps=1.e-6, maxit=200, weno=1, rotated_template=0, nsnx=10, nsnz=10, n_secondary=3, n_tertiary=3, radius_factor_tertiary=3.0, tt_from_rp=0, fsm_gpu=False) ->)

Parameters:

x (np.ndarray) – node coordinates along x, evenly spaced
z (np.ndarray) – node coordinates along z, evenly spaced
n_threads (int) – number of threads for raytracing (default is 1)
cell_slowness (bool) – slowness defined for cells (True) or nodes (False) (default is 1)
method (string) –
raytracing method (default is SPM)
- ’FSM’ : fast sweeping method
- ’SPM’ : shortest path method
- ’DSPM’ : dynamic shortest path method
dx and dz need not be equal, but the spacing along a given axis must be constant; a ValueError is raised otherwise.
aniso (string) –
type of anisotropy (implemented only for the SPM method)
- ’iso’ : isotropic medium
- ’elliptical’ : elliptical anisotropy
- ’tilted_elliptical’ : tilted elliptical anisotropy
- ’vti_psv’ : vertical transverse isotropy, P and SV waves
- ’vti_sh’ : vertical transverse isotropy, SH waves
- ’tti_psv’ : tilted transverse isotropy, P and SV waves
- ’tti_sh’ : tilted transverse isotropy, SH waves
- ’weakly_anelliptical’ : Weakly-Anelliptical formulation of B. Rommel
eps (double) – relative convergence criterion (FSM): the sweeps stop once the mean change in traveltime per node falls below this fraction of the traveltime range of the solution, so the same value behaves the same whatever units the model is expressed in (default is 1e-6)
maxit (int) – max number of sweeping iterations (FSM) (default is 200)
weno (bool) – use 3rd order weighted essentially non-oscillatory operator (FSM) (default is True)
rotated_template (bool) – use rotated templates (FSM, ignored when fsm_gpu is True)
nsnx (int) – number of secondary nodes in x (SPM) (default is 10)
nsnz (int) – number of secondary nodes in z (SPM) (default is 10)
n_secondary (int) – number of secondary nodes (DSPM) (default is 3)
n_tertiary (int) – number of tertiary nodes (DSPM) (default is 3)
radius_factor_tertiary (double) – multiplication factor used to compute radius of sphere around source that includes tertiary nodes (DSPM). The radius is the average edge length multiplied by this factor (default is 3)
tt_from_rp (bool) – compute traveltime using raypaths (available for FSM and DSPM only) (default is False)
fsm_gpu (bool) – use OpenCL implementation to run the Fast Sweeping Method on GPU (FSM only, default is False)

Notes

For raytracing in anisotropic media, the convention for inputting slowness depends on the model. For elliptical anisotropy, the method set_slowness is used to input horizontal slowness, while for weakly anelliptical anisotropy, the method is used to input vertical slowness.

compute_D(coord)

Return matrix of interpolation weights for velocity data points constraint

Parameters:: coord (np.ndarray with shape (npts, 2)) – coordinates of data points
Returns:: D – Matrix of interpolation weights
Return type:: scipy csr_array with shape (npts, nparams)

Note

In the current implementation, no check is made to see if the coordinates are on a node, edge, or corner.

compute_K(order=1)

Compute smoothing matrices

Parameters:: order (int) – order of smoothing operator, accept 1 or 2 (1 by default)
Returns:: Kx, Kz – matrices for derivatives along x & z
Return type:: tuple of csr_array

static data_kernel_straight_rays(Tx, Rx, grx, grz, aniso=False) → L

Raytracing with straight rays in 2D

Parameters:

Tx (np.ndarray) –
source coordinates, nTx by 2
- 1st column contains X coordinates,
- 2nd contains Z coordinates
Rx (np.ndarray) –
receiver coordinates, nTx by 2
- 1st column contains X coordinates,
- 2nd contains Z coordinates
grx (np.ndarray) – grid node coordinates along x
grz (np.ndarray) – grid node coordinates along z
aniso (bool) – compute L for elliptically anisotropic medium (True) or isotropic medium (False)

Returns:

L – data kernel matrix (tt = L @ slowness)

Return type:

scipy csr_array

Note

Tx and Rx should contain the same number of rows, each row corresponding to a source-receiver pair

dx

node separation along x

Type:: float

dz

node separation along x

Type:: float

get_grid_traveltimes(thread_no=0)

Obtain traveltimes computed at primary grid nodes

Parameters:: thread_no (int) – thread used to computed traveltimes (default is 0)
Returns:: tt
Return type:: np ndarray with shape (nx, nz)

get_niter()

Returns:: number of sweeping iterations performed by the last call to raytrace (FSM only, 0 for the other methods). When the WENO operator is used, this counts the first-order pass that precedes it, see get_niterw.
Return type:: int

Notes

A value equal to maxit means the sweeps ran out of iterations rather than reaching the convergence criterion, i.e. the traveltimes are not converged. A warning is then written to stderr by the solver.

When several sources are raytraced, the count is that of the source solved last, and with more than one thread it is whichever source finished last.

get_niterw()

Returns:: number of WENO sweeping iterations performed by the last call to raytrace (FSM with weno=1 only, 0 otherwise)
Return type:: int

Notes

The same caveats as for get_niter apply: a value equal to maxit means the WENO pass did not converge, and the count refers to the source solved last.

get_number_of_cells()

Returns:: number of cells in grid
Return type:: int

get_number_of_nodes()

Returns:: number of nodes in grid
Return type:: int

get_s0(hypo, slowness=None)

Return slowness at source points

Parameters:

hypo (np.ndarray with 5 columns) –
hypo holds source information, i.e.
- 1st column is event ID number
- 2nd column is origin time
- 3rd column is source easting (X)
- 4th column is source elevation (Z)
slowness (np ndarray with shape (nx, nz) (optional)) – slowness at grid nodes or cells (depending on cell_slowness) slowness may also have been flattened (with default ‘C’ order)

Returns:

s0 – slowness at source points

Return type:

np.ndarray

get_slowness()

Returns slowness of grid

Returns:: slowness
Return type:: np ndarray with shape (nx, nz)

Notes

Shape size will vary depending on slowness attribution to cells or nodes

is_outside(pts)

Check if points are outside grid

Parameters:: pts (np ndarray with shape (npts, 3)) – coordinates of points to check
Returns:: True if at least one point outside grid
Return type:: bool

n_threads

number of threads for raytracing

Type:: int

nparams

total number of parameters for grid

Type:: int

raytrace(source, rcv, slowness=None, xi=None, theta=None, Vp0=None, Vs0=None, delta=None, epsilon=None, gamma=None, thread_no=None, aggregate_src=False, compute_L=False, return_rays=False) → tt, rays, L

Perform raytracing

Parameters:

source (2D np.ndarray with 2 or 3 columns) – see notes below
rcv (2D np.ndarray with 2 columns) – Columns correspond to x, y and z coordinates
slowness (np ndarray with shape (nx, nz) (None by default)) – slowness at grid nodes or cells (depending on cell_slowness) slowness may also have been flattened (with default ‘C’ order) if None, slowness must have been assigned previously
xi (np ndarray with shape (nx, nz) (None by default)) – xi at grid cells (only for SPM & cell_slowness=True) xi may also have been flattened (with default ‘C’ order) if None, xi must have been assigned previously
theta (np ndarray with shape (nx, nz) (None by default)) – theta at grid cells (only for SPM & cell_slowness=True) theta may also have been flattened (with default ‘C’ order) if None, theta must have been assigned previously
Vp0 (np ndarray with shape (nx, nz) (None by default)) – Vp0 at grid cells (only for SPM & cell_slowness=True) Vp0 may also have been flattened (with default ‘C’ order) if None, Vp0 must have been assigned previously
Vs0 (np ndarray with shape (nx, nz) (None by default)) – Vs0 at grid cells (only for SPM & cell_slowness=True) Vs0 may also have been flattened (with default ‘C’ order) if None, Vs0 must have been assigned previously
delta (np ndarray with shape (nx, nz) (None by default)) – delta at grid cells (only for SPM & cell_slowness=True) delta may also have been flattened (with default ‘C’ order) if None, delta must have been assigned previously
epsilon (np ndarray with shape (nx, nz) (None by default)) – epsilon at grid cells (only for SPM & cell_slowness=True) epsilon may also have been flattened (with default ‘C’ order) if None, epsilon must have been assigned previously
gamma (np ndarray with shape (nx, nz) (None by default)) – gamma at grid cells (only for SPM & cell_slowness=True) gamma may also have been flattened (with default ‘C’ order) if None, gamma must have been assigned previously
thread_no (int (None by default)) – Perform calculations in thread number “thread_no” if None, attempt to run in parallel if warranted by number of sources and value of n_threads in constructor
aggregate_src (bool (False by default)) – if True, all source coordinates belong to a single event

compute_L (bool (False by default)) –

Compute matrices of partial derivative of travel time w/r to the medium parameters. L holds one block of ncells columns per parameter, in the order the setters take them:

aniso	blocks of columns
iso	slowness
elliptical	slowness, xi
vti_sh	Vs0, gamma
tilted_elliptical	slowness, xi, tilt angle
tti_sh	Vs0, gamma, tilt angle
weakly_anelliptical	slowness, s2, s4
vti_psv	Vp0, Vs0, epsilon, delta
tti_psv	Vp0, Vs0, epsilon, delta, tilt angle

return_rays (bool (False by default)) – Return raypaths

Returns:

tt (np.ndarray) – travel times for the appropriate source-rcv (see Notes below)
rays (list of np.ndarray) – Coordinates of segments forming raypaths (if return_rays is True)
L (scipy csr_array) – Matrix of partial derivative of travel time w/r to slowness

Notes

If source has 2 columns:

Columns correspond to x and z coordinates
Origin time (t0) is 0 for all points

If source has 3 columns:

1st column corresponds to origin times
2nd & 3rd columns correspond to x and z coordinates

source and rcv can contain the same number of rows, each row corresponding to a source-receiver pair, or the number of rows may differ if aggregate_src is True or if all rows in source are identical.

set_Vp0(v)

Assign vertical Vp to grid (VTI medium)

Parameters:: v (np ndarray with shape (nx, nz)) – v may also have been flattened (with default ‘C’ order)

set_Vs0(v)

Assign vertical Vs to grid (VTI medium)

Parameters:: v (np ndarray with shape (nx, nz)) – v may also have been flattened (with default ‘C’ order)

set_delta(d)

Assign Thomsen delta parameter to grid (VTI medium, P-SV waves)

Parameters:: d (np ndarray with shape (nx, nz)) – d may also have been flattened (with default ‘C’ order)

set_epsilon(e)

Assign Thomsen epsilon parameter to grid (VTI medium, P-SV waves)

Parameters:: e (np ndarray with shape (nx, nz)) – e may also have been flattened (with default ‘C’ order)

set_gamma(g)

Assign Thomsen gamma parameter to grid (VTI medium, SH waves)

Parameters:: g (np ndarray with shape (nx, nz)) – g may also have been flattened (with default ‘C’ order)

set_phase(phase)

Select the wave to model in a transversely isotropic medium

Parameters:: phase (str or int) – ‘qP’ for the quasi-compressional wave, ‘qSV’ for the quasi-shear one. The integers the C++ setPhase() takes are accepted as well, 1 for qP and anything else for qSV.

Notes

Only the ‘vti_psv’ and ‘tti_psv’ media describe both waves; the others raise. The qP wave is the one modelled until this is called.

set_s2(g)

Assign weakly anelliptical parameter s2

Parameters:: g (np ndarray with shape (nx, nz)) – g may also have been flattened (with default ‘C’ order)

set_s4(g)

Assign weakly anelliptical parameter s4

Parameters:: g (np ndarray with shape (nx, nz)) – g may also have been flattened (with default ‘C’ order)

set_slowness(slowness)

Assign slowness to grid

Parameters:: slowness (np ndarray with shape (nx, nz)) – slowness may also have been flattened (with default ‘C’ order)

set_tilt_angle(theta)

Assign anisotropy tilt angle to grid

Applies to the ‘tilted_elliptical’, ‘tti_psv’ and ‘tti_sh’ media. The angle is measured in radians; the symmetry axis lies at -theta from the vertical.

Parameters:: theta (np ndarray with shape (nx, nz)) – theta may also have been flattened (with default ‘C’ order)

set_traveltime_from_raypath(ttrp)

Set option to compute traveltime using raypath

Parameters:: ttrp (bool) – option value

set_use_thread_pool(use_thread_pool)

Set option to use thread pool instead of parallel loop

Parameters:: use_thread_pool (bool) – option value

set_velocity(velocity)

Assign velocity to grid

Parameters:: velocity (np ndarray with shape (nx, nz)) – velocity may also have been flattened (with default ‘C’ order)

set_xi(xi)

Assign elliptical anisotropy ratio to grid

Parameters:: xi (np ndarray with shape (nx, nz)) – xi may also have been flattened (with default ‘C’ order)

shape

number of parameters along each dimension

Type:: list of int

to_vtk(fields, filename)

Save grid variables and/or raypaths to VTK format

Parameters:

fields (dict) – dict of variables to save to file. Variables should be np.ndarray of size equal to either the number of nodes of the number of cells of the grid, or a list of raypath coordinates.
filename (str) – Name of file without extension for saving (extension vtr will be added). Raypaths are saved in separate files, and filename will be appended by the dict key and have a vtp extension.

Notes

VTK files can be visualized with Paraview (https://www.paraview.org)

x

node coordinates along x

Type:: np.ndarray

z

node coordinates along z

Type:: np.ndarray

class ttcrpy.rgrid.Grid2d_f

class to perform raytracing with 2D rectilinear grids (single precision)

Same as Grid2d_d but uses float32 for coordinates, slowness, traveltimes, and raypaths. Use the Grid2d factory to select precision via dtype=.

Constructor:

Grid2d_f(x, z, n_threads=1, cell_slowness=1, method=’SPM’, aniso=’iso’, eps=1.e-6, maxit=200, weno=1, rotated_template=0, nsnx=10, nsnz=10, n_secondary=3, n_tertiary=3, radius_factor_tertiary=3.0, tt_from_rp=0, fsm_gpu=False) -> Grid2d_f

Parameters:

x (np.ndarray (float32)) – node coordinates along x
z (np.ndarray (float32)) – node coordinates along z
Grid2d_d. (All other parameters are identical to)

compute_D(pts)

Return matrix of interpolation weights for velocity data points constraint

Parameters:: pts (np.ndarray with shape (npts, 2)) – coordinates of data points
Returns:: D – Matrix of interpolation weights
Return type:: scipy csr_array with shape (npts, nparams)

Note

In the current implementation, no check is made to see if the points are on a node, edge, or corner.

compute_K(order=1)

Compute smoothing matrices

Parameters:: order (int) – order of smoothing operator, accept 1 or 2 (1 by default)
Returns:: Kx, Kz – matrices for derivatives along x & z
Return type:: tuple of csr_array

dx

node separation along x

Type:: float

dz

node separation along x

Type:: float

get_grid_traveltimes(thread_no=0)

Obtain traveltimes computed at primary grid nodes

Parameters:: thread_no (int) – thread used to computed traveltimes (default is 0)
Returns:: tt
Return type:: np ndarray with shape (nx, nz)

get_niter()

Returns:: number of sweeping iterations performed by the last call to raytrace (FSM only, 0 for the other methods). When the WENO operator is used, this counts the first-order pass that precedes it, see get_niterw.
Return type:: int

Notes

A value equal to maxit means the sweeps ran out of iterations rather than reaching the convergence criterion, i.e. the traveltimes are not converged. A warning is then written to stderr by the solver.

When several sources are raytraced, the count is that of the source solved last, and with more than one thread it is whichever source finished last.

get_niterw()

Returns:: number of WENO sweeping iterations performed by the last call to raytrace (FSM with weno=1 only, 0 otherwise)
Return type:: int

Notes

The same caveats as for get_niter apply: a value equal to maxit means the WENO pass did not converge, and the count refers to the source solved last.

get_number_of_cells()

Returns:: number of cells in grid
Return type:: int

get_number_of_nodes()

Returns:: number of nodes in grid
Return type:: int

get_s0(hypo, slowness=None)

Return slowness at source points

Parameters:

hypo (np.ndarray with 5 columns) –
hypo holds source information, i.e.
- 1st column is event ID number
- 2nd column is origin time
- 3rd column is source easting (X)
- 4th column is source elevation (Z)
slowness (np ndarray with shape (nx, nz) (optional)) – slowness at grid nodes or cells (depending on cell_slowness) slowness may also have been flattened (with default ‘C’ order)

Returns:

s0 – slowness at source points

Return type:

np.ndarray

get_slowness()

Returns slowness of grid

Returns:: slowness
Return type:: np ndarray with shape (nx, nz)

Notes

Shape size will vary depending on slowness attribution to cells or nodes

is_outside(pts)

Check if points are outside grid

Parameters:: pts (np ndarray with shape (npts, 2)) – coordinates of points to check
Returns:: True if at least one point outside grid
Return type:: bool

n_threads

number of threads for raytracing

Type:: int

nparams

total number of parameters for grid

Type:: int

raytrace(source, rcv, slowness=None, xi=None, theta=None, Vp0=None, Vs0=None, delta=None, epsilon=None, gamma=None, thread_no=None, aggregate_src=False, compute_L=False, return_rays=False) → tt, rays, L

Perform raytracing

Parameters:

source (2D np.ndarray with 2 or 3 columns) – see notes below
rcv (2D np.ndarray with 2 columns) – Columns correspond to x and z coordinates
slowness (np ndarray with shape (nx, nz) (None by default)) – slowness at grid nodes or cells (depending on cell_slowness) slowness may also have been flattened (with default ‘C’ order) if None, slowness must have been assigned previously
xi (np ndarray with shape (nx, nz) (None by default)) – xi at grid cells (only for SPM & cell_slowness=True) xi may also have been flattened (with default ‘C’ order) if None, xi must have been assigned previously
theta (np ndarray with shape (nx, nz) (None by default)) – theta at grid cells (only for SPM & cell_slowness=True) theta may also have been flattened (with default ‘C’ order) if None, theta must have been assigned previously
Vp0 (np ndarray with shape (nx, nz) (None by default)) – Vp0 at grid cells (only for SPM & cell_slowness=True) Vp0 may also have been flattened (with default ‘C’ order) if None, Vp0 must have been assigned previously
Vs0 (np ndarray with shape (nx, nz) (None by default)) – Vs0 at grid cells (only for SPM & cell_slowness=True) Vs0 may also have been flattened (with default ‘C’ order) if None, Vs0 must have been assigned previously
delta (np ndarray with shape (nx, nz) (None by default)) – delta at grid cells (only for SPM & cell_slowness=True) delta may also have been flattened (with default ‘C’ order) if None, delta must have been assigned previously
epsilon (np ndarray with shape (nx, nz) (None by default)) – epsilon at grid cells (only for SPM & cell_slowness=True) epsilon may also have been flattened (with default ‘C’ order) if None, epsilon must have been assigned previously
gamma (np ndarray with shape (nx, nz) (None by default)) – gamma at grid cells (only for SPM & cell_slowness=True) gamma may also have been flattened (with default ‘C’ order) if None, gamma must have been assigned previously
thread_no (int (None by default)) – Perform calculations in thread number “thread_no” if None, attempt to run in parallel if warranted by number of sources and value of n_threads in constructor
aggregate_src (bool (False by default)) – if True, all source coordinates belong to a single event

compute_L (bool (False by default)) –

Compute matrices of partial derivative of travel time w/r to the medium parameters. L holds one block of ncells columns per parameter, in the order the setters take them:

aniso	blocks of columns
iso	slowness
elliptical	slowness, xi
vti_sh	Vs0, gamma
tilted_elliptical	slowness, xi, tilt angle
tti_sh	Vs0, gamma, tilt angle
weakly_anelliptical	slowness, s2, s4
vti_psv	Vp0, Vs0, epsilon, delta
tti_psv	Vp0, Vs0, epsilon, delta, tilt angle

return_rays (bool (False by default)) – Return raypaths

Returns:

tt (np.ndarray) – travel times for the appropriate source-rcv (see Notes below)
rays (list of np.ndarray) – Coordinates of segments forming raypaths (if return_rays is True)
L (scipy csr_array) – Matrix of partial derivative of travel time w/r to slowness

Notes

If source has 2 columns:

Columns correspond to x and z coordinates
Origin time (t0) is 0 for all points

If source has 3 columns:

1st column corresponds to origin times
2nd & 3rd columns correspond to x and z coordinates

source and rcv can contain the same number of rows, each row corresponding to a source-receiver pair, or the number of rows may differ if aggregate_src is True or if all rows in source are identical.

set_Vp0(v)

Assign vertical Vp to grid (VTI medium)

Parameters:: v (np ndarray with shape (nx, nz)) – v may also have been flattened (with default ‘C’ order)

set_Vs0(v)

Assign vertical Vs to grid (VTI medium)

Parameters:: v (np ndarray with shape (nx, nz)) – v may also have been flattened (with default ‘C’ order)

set_delta(d)

Assign Thomsen delta parameter to grid (VTI medium, P-SV waves)

Parameters:: d (np ndarray with shape (nx, nz)) – d may also have been flattened (with default ‘C’ order)

set_epsilon(e)

Assign Thomsen epsilon parameter to grid (VTI medium, P-SV waves)

Parameters:: e (np ndarray with shape (nx, nz)) – e may also have been flattened (with default ‘C’ order)

set_gamma(g)

Assign Thomsen gamma parameter to grid (VTI medium, SH waves)

Parameters:: g (np ndarray with shape (nx, nz)) – g may also have been flattened (with default ‘C’ order)

set_phase(phase)

Select the wave to model in a transversely isotropic medium

Parameters:: phase (str or int) – ‘qP’ for the quasi-compressional wave, ‘qSV’ for the quasi-shear one. The integers the C++ setPhase() takes are accepted as well, 1 for qP and anything else for qSV.

Notes

Only the ‘vti_psv’ and ‘tti_psv’ media describe both waves; the others raise. The qP wave is the one modelled until this is called.

set_s2(g)

Assign weakly anelliptical parameter s2

Parameters:: g (np ndarray with shape (nx, nz)) – g may also have been flattened (with default ‘C’ order)

set_s4(g)

Assign weakly anelliptical parameter s4

Parameters:: g (np ndarray with shape (nx, nz)) – g may also have been flattened (with default ‘C’ order)

set_slowness(slowness)

Assign slowness to grid

Parameters:: slowness (np ndarray with shape (nx, nz)) – slowness may also have been flattened (with default ‘C’ order)

set_tilt_angle(theta)

Assign anisotropy tilt angle to grid

Applies to the ‘tilted_elliptical’, ‘tti_psv’ and ‘tti_sh’ media. The angle is measured in radians; the symmetry axis lies at -theta from the vertical.

Parameters:: theta (np ndarray with shape (nx, nz)) – theta may also have been flattened (with default ‘C’ order)

set_traveltime_from_raypath(ttrp)

Set option to compute traveltime using raypath

Parameters:: ttrp (bool) – option value

set_use_thread_pool(use_thread_pool)

Set option to use thread pool instead of parallel loop

Parameters:: use_thread_pool (bool) – option value

set_velocity(velocity)

Assign velocity to grid

Parameters:: velocity (np ndarray with shape (nx, nz)) – velocity may also have been flattened (with default ‘C’ order)

set_xi(xi)

Assign elliptical anisotropy ratio to grid

Parameters:: xi (np ndarray with shape (nx, nz)) – xi may also have been flattened (with default ‘C’ order)

shape

number of parameters along each dimension

Type:: list of int

to_vtk(fields, filename)

Save grid variables and/or raypaths to VTK format

Parameters:

fields (dict) – dict of variables to save to file. Variables should be np.ndarray of size equal to either the number of nodes of the number of cells of the grid, or a list of raypath coordinates.
filename (str) – Name of file without extension for saving (extension vtr will be added). Raypaths are saved in separate files, and filename will be appended by the dict key and have a vtp extension.

Notes

VTK files can be visualized with Paraview (https://www.paraview.org)

x

node coordinates along x

Type:: np.ndarray

z

node coordinates along z

Type:: np.ndarray

ttcrpy.rgrid.Grid3d(x, y, z, ..., aniso='iso', ..., dtype=np.float64) → Grid3d_d or Grid3d_f

Factory that returns a 3D rectilinear-grid raytracer with the requested numeric precision.

Parameters:

x (array-like) – Node coordinates. Converted to the requested dtype automatically.
y (array-like) – Node coordinates. Converted to the requested dtype automatically.
z (array-like) – Node coordinates. Converted to the requested dtype automatically.
dtype (numpy dtype) – np.float64 (default) for double precision (Grid3d_d), np.float32 for single precision (Grid3d_f).
class. (All other parameters are forwarded unchanged to the underlying)

Return type:

Grid3d_d or Grid3d_f

class ttcrpy.rgrid.Grid3d_d

class to perform raytracing with 3D rectilinear grids (double precision)

Variables:

x (np.ndarray) – node coordinates along x
y (np.ndarray) – node coordinates along y
z (np.ndarray) – node coordinates along z
dx (float) – node separation along x
dy (float) – node separation along y
dz (float) – node separation along z
shape ((int, int, int)) – number of parameters along each dimension
nparams (int) – total number of parameters for grid
n_threads (int) – number of threads for raytracing
Constructor

Grid3d_d (Grid3d_d(x, y, z, n_threads=1, cell_slowness=1, method='FSM', tt_from_rp=1, interp_vel=0, eps=1.e-6, maxit=200, weno=1, nsnx=5, nsny=5, nsnz=5, n_secondary=2, n_tertiary=2, radius_factor_tertiary=3.0, translate_grid=False) ->) –

param x:

node coordinates along x, evenly spaced

type x:

np.ndarray

param y:

node coordinates along y, evenly spaced

type y:

np.ndarray

param z:

node coordinates along z, evenly spaced

type z:

np.ndarray

param n_threads:

number of threads for raytracing (default is 1)

type n_threads:

int

param cell_slowness:

slowness defined for cells (True) or nodes (False) (default is 1)

type cell_slowness:

bool

param method:

raytracing method (default is FSM)

’FSM’ : fast sweeping method

’SPM’ : shortest path method

’DSPM’ : dynamic shortest path

The three node spacings need not be equal for any of the methods. FSM picks a general per-axis stencil when they differ and a cheaper equal-spacing one when they do not, for the WENO stencil as well as the first-order one. The spacing along a given axis must be constant, however; a ValueError is raised otherwise.

type method:

string

param aniso:

type of anisotropy (SPM method and cell_slowness only)

’iso’ : isotropic medium
’elliptical’ : ellipsoidal anisotropy, axes aligned with the grid; set_slowness takes the vertical slowness and the two ratios are given with set_chi and set_psi
’vti_psv’ : vertical transverse isotropy, P and SV waves
’vti_sh’ : vertical transverse isotropy, SH waves
’weakly_anelliptical’ : Weakly-Anelliptical formulation of B. Rommel; set_slowness takes the vertical slowness

The tilted models of Grid2d have no 3D counterpart yet. The parameters of each model, and the order of the blocks of columns compute_L returns, are

aniso	setters, in the order the blocks appear
’iso’	set_slowness
’elliptical’	set_slowness, set_chi, set_psi
’vti_sh’	set_Vs0, set_gamma
’weakly_anelliptical’	set_slowness, set_s2, set_s4
’vti_psv’	set_Vp0, set_Vs0, set_epsilon, set_delta

type aniso:

string

param tt_from_rp:

compute traveltimes from raypaths (FSM or DSPM only) (default is 1)

type tt_from_rp:

bool

param interp_vel:

interpolate velocity instead of slowness at nodes (for cell_slowness == False or FSM) (defauls is False)

type interp_vel:

bool

param eps:

relative convergence criterion (FSM): the sweeps stop once the mean change in traveltime per node falls below this fraction of the traveltime range of the solution, so the same value behaves the same whatever units the model is expressed in (default is 1e-6)

type eps:

double

param maxit:

max number of sweeping iterations (FSM) (default is 200)

type maxit:

int

param weno:

use 3rd order weighted essentially non-oscillatory operator (FSM) (default is True)

type weno:

bool

param nsnx:

number of secondary nodes in x (SPM) (default is 5)

type nsnx:

int

param nsny:

number of secondary nodes in y (SPM) (default is 5)

type nsny:

int

param nsnz:

number of secondary nodes in z (SPM) (default is 5)

type nsnz:

int

param n_secondary:

number of secondary nodes (DSPM) (default is 2)

type n_secondary:

int

param n_tertiary:

number of tertiary nodes (DSPM) (default is 2)

type n_tertiary:

int

param radius_factor_tertiary:

multiplication factor used to compute radius of sphere around source that includes tertiary nodes (DSPM). The radius is the average edge length multiplied by this factor (default is 3)

type radius_factor_tertiary:

double

param translate_grid:

Translate the grid such that origin is (0, 0, 0) to perform computations, which may increase accuracy when large values, e.g. UTM coordinates, are used. When raytracing, src and rcv should be given in the original system, and output raypath coordinates are also given in the original system (default if False)

type translate_grid:

bool

param fsm_gpu:

Use OpenCL implementation to run the Fast Sweeping Method on GPU (default is False)

type fsm_gpu:

bool

static builder(filename, n_threads=1, method='FSM', tt_from_rp=1, interp_vel=0, eps=1.e-6, maxit=200, weno=1, nsnx=5, nsny=5, nsnz=5, n_secondary=2, n_tertiary=2, radius_factor_tertiary=3.0, translate_grid=0)

Build instance of Grid3d from VTK file

Parameters:

filename (str) – Name of file holding a vtkRectilinearGrid. The grid must have point or cell attribute named either ‘Slowness’, ‘slowness’, ‘Velocity’, ‘velocity’, or ‘P-wave velocity’
Constructor (Other parameters are defined in)

Returns:

grid – grid instance

Return type:

Grid3d

compute_D(coord)

Return matrix of interpolation weights for velocity data points constraint

Parameters:: coord (np.ndarray with shape (npts, 3)) – coordinates of data points
Returns:: D – Matrix of interpolation weights
Return type:: scipy csr_array with shape (npts, nparams)

Note

In the current implementation, no check is made to see if the coordinates are on a node, edge, face, or corner.

compute_H(source, rcv, slowness=None, full=True, radius_factor=4.0, thread_no=None) → tt, H

Traveltimes and the hypocentre-location Jacobian

H holds the partial derivatives of the arrival time with respect to the hypocentre parameters, one row per receiver. Raytracing is performed internally, so the traveltimes are returned along with H and a separate call to raytrace is not needed.

Parameters:

source (2D np.ndarray with 3, 4 or 5 columns) – see notes of raytrace
rcv (2D np.ndarray with 3 columns) – Columns correspond to x, y and z coordinates
slowness (np ndarray with shape (nx, ny, nz) (None by default)) – slowness at grid nodes or cells (depending on cell_slowness) if None, slowness must have been assigned previously
full (bool (True by default)) –
if True, H has four columns

[1, dT/dx, dT/dy, dT/dz]

the leading 1 being the derivative with respect to origin time. If False, H has the two columns [dT/dx, dT/dy].
radius_factor (double (4.0 by default)) – the take-off direction is measured where the walk back from the receiver first comes within radius_factor average edge lengths of the source. Closer than that, the traveltime field is radially degenerate about the source; much further out, the chord departs from the ray tangent. The error has a minimum in between, shallow on rectilinear grids and more pronounced on coarse meshes.
thread_no (int (None by default)) – thread number to use (a single source is then expected)

Returns:

tt (np.ndarray with shape (nrcv,)) – traveltimes
H (np.ndarray with shape (nrcv, 4) or (nrcv, 2)) – Jacobian

Notes

The spatial derivatives follow from dT/dx_s = -s(x_s) * e, with e the unit take-off direction at the source. e is obtained by descending the traveltime field from the receiver rather than from the raypath: the raypath endpoint convention differs between solvers, and its final segment is a noisy estimate of the tangent.

compute_K()

Compute smoothing matrices (2nd order derivative)

Returns:: Kx, Ky, Kz – matrices for derivatives along x, y, & z
Return type:: tuple of csr_array

static data_kernel_straight_rays(Tx, Rx, grx, gry, grz, centers) -> L, (xc, yc, zc)

Raytracing with straight rays in 3D

Parameters:

Tx (np.ndarray) –
source coordinates, nTx by 3
- 1st column contains X coordinates,
- 2nd contains Y coordinates
- 3rd contains Z coordinates
Rx (np.ndarray) –
receiver coordinates, nTx by 3
- 1st column contains X coordinates,
- 2nd contains Y coordinates
- 3rd contains Z coordinates
grx (np.ndarray) – grid node coordinates along x
gry (np.ndarray) – grid node coordinates along y
grz (np.ndarray) – grid node coordinates along z
centers (bool) – return coordinates of center of cells (False by default)

Returns:

L (scipy csr_array) – data kernel matrix (tt = L @ slowness)
(xc, yc, zc) (tuple of np.ndarray) – vectors of coordinates of center of cells

Note

Tx and Rx should contain the same number of rows, each row corresponding to a source-receiver pair

dx

node separation along x

Type:: float

dy

node separation along y

Type:: float

dz

node separation along z

Type:: float

get_grid_traveltimes(thread_no=0)

Obtain traveltimes computed at primary grid nodes

Parameters:: thread_no (int) – thread used to computed traveltimes (default is 0)
Returns:: tt – traveltimes
Return type:: np ndarray with shape (nx, ny, nz)

get_niter()

Returns:: number of sweeping iterations performed by the last call to raytrace (FSM only, 0 for the other methods). When the WENO operator is used, this counts the first-order pass that precedes it, see get_niterw.
Return type:: int

Notes

A value equal to maxit means the sweeps ran out of iterations rather than reaching the convergence criterion, i.e. the traveltimes are not converged. A warning is then written to stderr by the solver.

When several sources are raytraced, the count is that of the source solved last, and with more than one thread it is whichever source finished last.

get_niterw()

Returns:: number of WENO sweeping iterations performed by the last call to raytrace (FSM with weno=1 only, 0 otherwise)
Return type:: int

Notes

The same caveats as for get_niter apply: a value equal to maxit means the WENO pass did not converge, and the count refers to the source solved last.

get_number_of_cells()

Returns:: number of cells in grid
Return type:: int

get_number_of_nodes()

Returns:: number of nodes in grid
Return type:: int

get_s0(hypo, slowness=None)

Return slowness at source points

Parameters:

hypo (np.ndarray with 5 columns) –
hypo holds source information, i.e.
- 1st column is event ID number
- 2nd column is origin time
- 3rd column is source easting
- 4th column is source northing
- 5th column is source elevation
slowness (np ndarray with shape (nx, ny, nz) (optional)) – slowness at grid nodes or cells (depending on cell_slowness) slowness may also have been flattened (with default ‘C’ order)

Returns:

s0 – slowness at source points

Return type:

np.ndarray

get_slowness()

Returns slowness of grid

Returns:: slowness
Return type:: np ndarray with shape (nx, ny, nz)

Notes

Shape size will vary depending on slowness attribution to cells or nodes

ind(i, j, k)

Return node index

Parameters:

i (int) – index of node along x
j (int) – index of node along y
k (int) – index of node along z

Returns:

node index for a “flattened” grid

Return type:

int

indc(i, j, k)

return cell index

Parameters:

i (int) – index of cell along x
j (int) – index of cell along y
k (int) – index of cell along z

Returns:

cell index for a “flattened” grid

Return type:

int

is_outside(pts)

Check if points are outside grid

Parameters:: pts (np ndarray with shape (npts, 3)) – coordinates of points to check
Returns:: True if at least one point outside grid
Return type:: bool

is_using_gpu

whether the solve actually runs on the GPU

fsm_gpu asks for the OpenCL solvers; it does not guarantee them. The request is refused when no device is available, when initialisation fails, or when the grid is double precision and the device reports no cl_khr_fp64 – which is every Apple GPU, so np.float32 is what reaches the GPU there. A refusal falls back to the CPU and leaves the results correct, so this property is the only way to tell the two apart.

Type:: bool

n_threads

number of threads for raytracing

Type:: int

nparams

total number of parameters for grid

Type:: int

raytrace(source, rcv, slowness=None, thread_no=None, aggregate_src=False, compute_L=False, compute_M=False, return_rays=False)

raytrace(source, rcv, slowness=None, thread_no=None,: aggregate_src=False, compute_L=False, compute_M=False, return_rays=False) -> tt, rays, M, L

Perform raytracing

Parameters:

source (2D np.ndarray with 3, 4 or 5 columns) – see notes below
rcv (2D np.ndarray with 3 columns) – Columns correspond to x, y and z coordinates
slowness (np ndarray with shape (nx, ny, nz) (None by default)) – slowness at grid nodes or cells (depending on cell_slowness) slowness may also have been flattened (with default ‘C’ order) if None, slowness must have been assigned previously
thread_no (int (None by default)) – Perform calculations in thread number “thread_no” if None, attempt to run in parallel if warranted by number of sources and value of n_threads in constructor
aggregate_src (bool (False by default)) – if True, all source coordinates belong to a single event
compute_L (bool (False by default)) –
Compute matrices of partial derivative of travel time w/r to slowness. Requires slowness defined at cells; available for the FSM, SPM and DSPM.

L holds the length the ray spends in each cell, so L @ s is the traveltime the raypath integrates. For the SPM and DSPM that is the traveltime returned in tt, to machine precision. The FSM solves on the nodes, over the cell slownesses averaged onto them (see Grid3Drcfs), and integrates that interpolated field along the path, so its tt and L @ s differ by a discretization term – a few parts in 1e3 on a coarse grid, falling with the cell size.
compute_M (bool (False by default)) – Compute matrices of partial derivative of travel time w/r to velocity Note : compute_M and compute_L are mutually exclusive
return_rays (bool (False by default)) –
Return raypaths

Keep sources at least one cell away from the edges of the model. The FSM and DSPM rebuild raypaths by descending the traveltime gradient and stop once within one voxel diagonal of the source; a path that reaches a face of the model before then has no step left and raises ‘Error while computing raypaths: going outside grid’. Where the velocity increases towards a boundary the fastest path rides it, so a source near that boundary is approached along the face. The SPM follows the node parents recorded during the sweep and is not affected.

Returns:

tt (np.ndarray) – travel times for the appropriate source-rcv (see Notes below)
rays (list of np.ndarray) – Coordinates of segments forming raypaths (if return_rays is True)
M (list of csr_array) – matrices of partial derivative of travel time w/r to velocity. the number of matrices is equal to the number of sources
L (scipy csr_array) – Matrix of partial derivative of travel time w/r to slowness. if input argument source has 5 columns, L is a list of matrices and the number of matrices is equal to the number of sources otherwise, L is a single csr_array

Notes

If source has 3 columns:

Columns correspond to x, y and z coordinates
Origin time (t0) is 0 for all points

If source has 4 columns:

1st column corresponds to origin times
2nd, 3rd & 4th columns correspond to x, y and z coordinates

If source has 5 columns:

1st column corresponds to event ID
2nd column corresponds to origin times
3rd, 4th & 5th columns correspond to x, y and z coordinates

For the latter case (5 columns), source and rcv should contain the same number of rows, each row corresponding to a source-receiver pair. For the 2 other cases, source and rcv can contain the same number of rows, each row corresponding to a source-receiver pair, or the number of rows may differ if aggregate_src is True or if all rows in source are identical.

set_Vp0(v)

Assign vertical P-wave velocity (transversely isotropic medium) to grid

Parameters:: v (np ndarray with shape (nx, ny, nz)) – v may also have been flattened (with default ‘C’ order)

set_Vs0(v)

Assign vertical S-wave velocity (transversely isotropic medium) to grid

Parameters:: v (np ndarray with shape (nx, ny, nz)) – v may also have been flattened (with default ‘C’ order)

set_chi(chi)

Assign elliptical anisotropy ratio \(\chi = s_x/s_z\) to grid

Parameters:: chi (np ndarray with shape (nx, ny, nz)) – chi may also have been flattened (with default ‘C’ order)

set_delta(delta)

Assign Thomsen’s parameter \(\delta\) to grid

Parameters:: delta (np ndarray with shape (nx, ny, nz)) – delta may also have been flattened (with default ‘C’ order)

set_epsilon(epsilon)

Assign Thomsen’s parameter \(\epsilon\) to grid

Parameters:: epsilon (np ndarray with shape (nx, ny, nz)) – epsilon may also have been flattened (with default ‘C’ order)

set_gamma(gamma)

Assign Thomsen’s parameter \(\gamma\) to grid

Parameters:: gamma (np ndarray with shape (nx, ny, nz)) – gamma may also have been flattened (with default ‘C’ order)

set_phase(phase)

Select the wave to model in a transversely isotropic medium

Parameters:: phase (str or int) – ‘qP’ for the quasi-compressional wave, ‘qSV’ for the quasi-shear one. The integers the C++ setPhase() takes are accepted as well, 1 for qP and anything else for qSV.

Notes

Only the ‘vti_psv’ medium describes both waves; the others raise. The qP wave is the one modelled until this is called.

set_psi(psi)

Assign elliptical anisotropy ratio \(\psi = s_y/s_z\) to grid

Parameters:: psi (np ndarray with shape (nx, ny, nz)) – psi may also have been flattened (with default ‘C’ order)

set_s2(s2)

Assign second-order anisotropy coefficient (weakly anelliptical medium) to grid

Parameters:: s2 (np ndarray with shape (nx, ny, nz)) – s2 may also have been flattened (with default ‘C’ order)

set_s4(s4)

Assign fourth-order anisotropy coefficient (weakly anelliptical medium) to grid

Parameters:: s4 (np ndarray with shape (nx, ny, nz)) – s4 may also have been flattened (with default ‘C’ order)

set_slowness(slowness)

Assign slowness to grid

Parameters:: slowness (np ndarray with shape (nx, ny, nz)) – slowness may also have been flattened (with default ‘C’ order)

set_traveltime_from_raypath(ttrp)

Set option to compute traveltime using raypath

Parameters:: ttrp (bool) – option value

set_use_thread_pool(use_thread_pool)

Set option to use thread pool instead of parallel loop

Parameters:: use_thread_pool (bool) – option value

set_velocity(velocity)

Assign velocity to grid

Parameters:: velocity (np ndarray with shape (nx, ny, nz)) – velocity may also have been flattened (with default ‘C’ order)

shape

number of parameters along each dimension

Type:: list of int

to_vtk(fields, filename)

Save grid variables and/or raypaths to VTK format

Parameters:

fields (dict) – dict of variables to save to file. Variables should be np.ndarray of size equal to either the number of nodes of the number of cells of the grid, or a list of raypath coordinates.
filename (str) – Name of file without extension for saving (extension vtr will be added). Raypaths are saved in separate files, and filename will be appended by the dict key and have a vtp extension.

Notes

VTK files can be visualized with Paraview (https://www.paraview.org)

x

node coordinates along x

Type:: np.ndarray

y

node coordinates along y

Type:: np.ndarray

z

node coordinates along z

Type:: np.ndarray

class ttcrpy.rgrid.Grid3d_f

class to perform raytracing with 3D rectilinear grids (single precision)

Same as Grid3d_d but uses float32 for coordinates, slowness, traveltimes, and raypaths. Use the Grid3d factory to select precision via dtype=.

Constructor:

Grid3d_f(x, y, z, n_threads=1, cell_slowness=1, method=’FSM’, aniso=’iso’, tt_from_rp=1, interp_vel=0, eps=1.e-6, maxit=200, weno=1, nsnx=5, nsny=5, nsnz=5, n_secondary=2, n_tertiary=2, radius_factor_tertiary=3.0, translate_grid=False, fsm_gpu=False) -> Grid3d_f

Parameters:

x (np.ndarray (float32)) – node coordinates along x
y (np.ndarray (float32)) – node coordinates along y
z (np.ndarray (float32)) – node coordinates along z
Grid3d_d. (All other parameters are identical to)

static builder(filename, ...) → Grid3d_f

Build instance of Grid3d_f from VTK file (single precision).

Parameters:

filename (str) – Name of file holding a vtkRectilinearGrid.
Constructor. (Other parameters are defined in the)

compute_D(pts)

Return matrix of interpolation weights for velocity data points constraint

Parameters:: pts (np.ndarray with shape (npts, 3)) – coordinates of data points
Returns:: D – Matrix of interpolation weights
Return type:: scipy csr_array with shape (npts, nparams)

Note

In the current implementation, no check is made to see if the points are on a node, edge, or corner.

compute_H(source, rcv, slowness=None, full=True, radius_factor=4.0, thread_no=None) → tt, H

Traveltimes and the hypocentre-location Jacobian

H holds the partial derivatives of the arrival time with respect to the hypocentre parameters, one row per receiver. Raytracing is performed internally, so the traveltimes are returned along with H and a separate call to raytrace is not needed.

Parameters:

source (2D np.ndarray with 3, 4 or 5 columns) – see notes of raytrace
rcv (2D np.ndarray with 3 columns) – Columns correspond to x, y and z coordinates
slowness (np ndarray with shape (nx, ny, nz) (None by default)) – slowness at grid nodes or cells (depending on cell_slowness) if None, slowness must have been assigned previously
full (bool (True by default)) –
if True, H has four columns

[1, dT/dx, dT/dy, dT/dz]

the leading 1 being the derivative with respect to origin time. If False, H has the two columns [dT/dx, dT/dy].
radius_factor (double (4.0 by default)) – the take-off direction is measured where the walk back from the receiver first comes within radius_factor average edge lengths of the source. Closer than that, the traveltime field is radially degenerate about the source; much further out, the chord departs from the ray tangent. The error has a minimum in between, shallow on rectilinear grids and more pronounced on coarse meshes.
thread_no (int (None by default)) – thread number to use (a single source is then expected)

Returns:

tt (np.ndarray with shape (nrcv,)) – traveltimes
H (np.ndarray with shape (nrcv, 4) or (nrcv, 2)) – Jacobian

Notes

The spatial derivatives follow from dT/dx_s = -s(x_s) * e, with e the unit take-off direction at the source. e is obtained by descending the traveltime field from the receiver rather than from the raypath: the raypath endpoint convention differs between solvers, and its final segment is a noisy estimate of the tangent.

compute_K()

Compute smoothing matrices (2nd order derivative)

Returns:: Kx, Ky, Kz – matrices for derivatives along x, y, & z
Return type:: tuple of csr_array

dx

node separation along x

Type:: float

dy

node separation along y

Type:: float

dz

node separation along z

Type:: float

get_grid_traveltimes(thread_no=0)

Obtain traveltimes computed at primary grid nodes

Parameters:: thread_no (int) – thread used to compute traveltimes (default is 0)
Returns:: tt – traveltimes
Return type:: np ndarray with shape (nx, ny, nz)

get_niter()

Returns:: number of sweeping iterations performed by the last call to raytrace (FSM only, 0 for the other methods). When the WENO operator is used, this counts the first-order pass that precedes it, see get_niterw.
Return type:: int

Notes

A value equal to maxit means the sweeps ran out of iterations rather than reaching the convergence criterion, i.e. the traveltimes are not converged. A warning is then written to stderr by the solver.

When several sources are raytraced, the count is that of the source solved last, and with more than one thread it is whichever source finished last.

get_niterw()

Returns:: number of WENO sweeping iterations performed by the last call to raytrace (FSM with weno=1 only, 0 otherwise)
Return type:: int

Notes

The same caveats as for get_niter apply: a value equal to maxit means the WENO pass did not converge, and the count refers to the source solved last.

get_number_of_cells()

Returns:: number of cells in grid
Return type:: int

get_number_of_nodes()

Returns:: number of nodes in grid
Return type:: int

get_s0(hypo, slowness=None)

Return slowness at source points

Parameters:

hypo (np.ndarray with 5 columns) –
hypo holds source information, i.e.
- 1st column is event ID number
- 2nd column is origin time
- 3rd column is source easting
- 4th column is source northing
- 5th column is source elevation
slowness (np ndarray with shape (nx, ny, nz) (optional)) – slowness at grid nodes or cells (depending on cell_slowness) slowness may also have been flattened (with default ‘C’ order)

Returns:

s0 – slowness at source points

Return type:

np.ndarray

get_slowness()

Returns slowness of grid

Returns:: slowness
Return type:: np ndarray with shape (nx, ny, nz)

Notes

Shape size will vary depending on slowness attribution to cells or nodes

ind(i, j, k)

Return node index

Parameters:

i (int) – index of node along x
j (int) – index of node along y
k (int) – index of node along z

Returns:

node index for a “flattened” grid

Return type:

int

indc(i, j, k)

return cell index

Parameters:

i (int) – index of cell along x
j (int) – index of cell along y
k (int) – index of cell along z

Returns:

cell index for a “flattened” grid

Return type:

int

is_outside(pts)

Check if points are outside grid

Parameters:: pts (np ndarray with shape (npts, 3)) – coordinates of points to check
Returns:: True if at least one point outside grid
Return type:: bool

is_using_gpu

whether the solve actually runs on the GPU

fsm_gpu asks for the OpenCL solvers; it does not guarantee them. The request is refused when no device is available, when initialisation fails, or when the grid is double precision and the device reports no cl_khr_fp64 – which is every Apple GPU, so np.float32 is what reaches the GPU there. A refusal falls back to the CPU and leaves the results correct, so this property is the only way to tell the two apart.

Type:: bool

n_threads

number of threads for raytracing

Type:: int

nparams

total number of parameters for grid

Type:: int

raytrace(source, rcv, slowness=None, thread_no=None, aggregate_src=False, compute_L=False, compute_M=False, return_rays=False)

raytrace(source, rcv, slowness=None, thread_no=None,: aggregate_src=False, compute_L=False, compute_M=False, return_rays=False) -> tt, rays, M, L

Perform raytracing

Parameters:

source (2D np.ndarray with 3, 4 or 5 columns) – see notes below
rcv (2D np.ndarray with 3 columns) – Columns correspond to x, y and z coordinates
slowness (np ndarray with shape (nx, ny, nz) (None by default)) – slowness at grid nodes or cells (depending on cell_slowness) slowness may also have been flattened (with default ‘C’ order) if None, slowness must have been assigned previously
thread_no (int (None by default)) – Perform calculations in thread number “thread_no” if None, attempt to run in parallel if warranted by number of sources and value of n_threads in constructor
aggregate_src (bool (False by default)) – if True, all source coordinates belong to a single event
compute_L (bool (False by default)) –
Compute matrices of partial derivative of travel time w/r to slowness. Requires slowness defined at cells; available for the FSM, SPM and DSPM.

L holds the length the ray spends in each cell, so L @ s is the traveltime the raypath integrates. For the SPM and DSPM that is the traveltime returned in tt, to machine precision. The FSM solves on the nodes, over the cell slownesses averaged onto them (see Grid3Drcfs), and integrates that interpolated field along the path, so its tt and L @ s differ by a discretization term – a few parts in 1e3 on a coarse grid, falling with the cell size.
compute_M (bool (False by default)) – Compute matrices of partial derivative of travel time w/r to velocity Note : compute_M and compute_L are mutually exclusive
return_rays (bool (False by default)) –
Return raypaths

Keep sources at least one cell away from the edges of the model. The FSM and DSPM rebuild raypaths by descending the traveltime gradient and stop once within one voxel diagonal of the source; a path that reaches a face of the model before then has no step left and raises ‘Error while computing raypaths: going outside grid’. Where the velocity increases towards a boundary the fastest path rides it, so a source near that boundary is approached along the face. The SPM follows the node parents recorded during the sweep and is not affected.

Returns:

tt (np.ndarray) – travel times for the appropriate source-rcv (see Notes below)
rays (list of np.ndarray) – Coordinates of segments forming raypaths (if return_rays is True)
M (list of csr_array) – matrices of partial derivative of travel time w/r to velocity. the number of matrices is equal to the number of sources
L (scipy csr_array) – Matrix of partial derivative of travel time w/r to slowness. if input argument source has 5 columns, L is a list of matrices and the number of matrices is equal to the number of sources otherwise, L is a single csr_array

Notes

If source has 3 columns:

Columns correspond to x, y and z coordinates
Origin time (t0) is 0 for all points

If source has 4 columns:

1st column corresponds to origin times
2nd, 3rd & 4th columns correspond to x, y and z coordinates

If source has 5 columns:

1st column corresponds to event ID
2nd column corresponds to origin times
3rd, 4th & 5th columns correspond to x, y and z coordinates

For the latter case (5 columns), source and rcv should contain the same number of rows, each row corresponding to a source-receiver pair. For the 2 other cases, source and rcv can contain the same number of rows, each row corresponding to a source-receiver pair, or the number of rows may differ if aggregate_src is True or if all rows in source are identical.

set_Vp0(v)

Assign vertical P-wave velocity (transversely isotropic medium) to grid

Parameters:: v (np ndarray with shape (nx, ny, nz)) – v may also have been flattened (with default ‘C’ order)

set_Vs0(v)

Assign vertical S-wave velocity (transversely isotropic medium) to grid

Parameters:: v (np ndarray with shape (nx, ny, nz)) – v may also have been flattened (with default ‘C’ order)

set_chi(chi)

Assign elliptical anisotropy ratio \(\chi = s_x/s_z\) to grid

Parameters:: chi (np ndarray with shape (nx, ny, nz)) – chi may also have been flattened (with default ‘C’ order)

set_delta(delta)

Assign Thomsen’s parameter \(\delta\) to grid

Parameters:: delta (np ndarray with shape (nx, ny, nz)) – delta may also have been flattened (with default ‘C’ order)

set_epsilon(epsilon)

Assign Thomsen’s parameter \(\epsilon\) to grid

Parameters:: epsilon (np ndarray with shape (nx, ny, nz)) – epsilon may also have been flattened (with default ‘C’ order)

set_gamma(gamma)

Assign Thomsen’s parameter \(\gamma\) to grid

Parameters:: gamma (np ndarray with shape (nx, ny, nz)) – gamma may also have been flattened (with default ‘C’ order)

set_phase(phase)

Select the wave to model in a transversely isotropic medium

Parameters:: phase (str or int) – ‘qP’ for the quasi-compressional wave, ‘qSV’ for the quasi-shear one. The integers the C++ setPhase() takes are accepted as well, 1 for qP and anything else for qSV.

Notes

Only the ‘vti_psv’ medium describes both waves; the others raise. The qP wave is the one modelled until this is called.

set_psi(psi)

Assign elliptical anisotropy ratio \(\psi = s_y/s_z\) to grid

Parameters:: psi (np ndarray with shape (nx, ny, nz)) – psi may also have been flattened (with default ‘C’ order)

set_s2(s2)

Assign second-order anisotropy coefficient (weakly anelliptical medium) to grid

Parameters:: s2 (np ndarray with shape (nx, ny, nz)) – s2 may also have been flattened (with default ‘C’ order)

set_s4(s4)

Assign fourth-order anisotropy coefficient (weakly anelliptical medium) to grid

Parameters:: s4 (np ndarray with shape (nx, ny, nz)) – s4 may also have been flattened (with default ‘C’ order)

set_slowness(slowness)

Assign slowness to grid

Parameters:: slowness (np ndarray with shape (nx, ny, nz)) – slowness may also have been flattened (with default ‘C’ order)

set_traveltime_from_raypath(ttrp)

Set option to compute traveltime using raypath

Parameters:: ttrp (bool) – option value

set_use_thread_pool(use_thread_pool)

Set option to use thread pool instead of parallel loop

Parameters:: use_thread_pool (bool) – option value

set_velocity(velocity)

Assign velocity to grid

Parameters:: velocity (np ndarray with shape (nx, ny, nz)) – velocity may also have been flattened (with default ‘C’ order)

shape

number of parameters along each dimension

Type:: list of int

to_vtk(fields, filename)

Save grid variables and/or raypaths to VTK format

Parameters:

fields (dict) – dict of variables to save to file. Variables should be np.ndarray of size equal to either the number of nodes or the number of cells of the grid, or a list of raypath coordinates.
filename (str) – Name of file without extension for saving (extension vtr will be added). Raypaths are saved in separate files, and filename will be appended by the dict key and have a vtp extension.

Notes

VTK files can be visualized with Paraview (https://www.paraview.org)

x

node coordinates along x

Type:: np.ndarray

y

node coordinates along y

Type:: np.ndarray

z

node coordinates along z

Type:: np.ndarray

ttcrpy.rgrid.set_verbose(v)

Set verbosity level for C++ code

Parameters:: v (int) – verbosity level