Module tmesh

Raytracing on unstructured triangular and tetrahedral meshes

This module contains two classes to perform traveltime computation and raytracing on unstructured meshes:

• Mesh2d for 2D media

• Mesh3d 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:

1. slowness constant within the voxels of the mesh (the default)

2. slowness defined at nodes of the mesh

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

class ttcrpy.tmesh.Mesh2d

class to perform raytracing with triangular meshes

Variables:

• nparams (int) – total number of parameters for grid

• n_threads (int) – number of threads for raytracing

• Constructor

• Mesh2d (Mesh2d(nodes, triangles, n_threads=1, cell_slowness=1, method='FSM', aniso='iso', eps=1e-6, maxit=200, process_obtuse=1, n_secondary=5, n_tertiary=2, radius_factor_tertiary=2, tt_from_rp=0) ->)

Parameters:

• nodes (np.ndarray with shape (nnodes, 2)) – node coordinates

• triangles (np.ndarray of int with shape (ntriangles, 3)) – indices of nodes forming the triangles

• 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 FSM)

  • ’FSM’ : fast sweeping method

  • ’SPM’ : shortest path method

  • ’DSPM’ : dynamic shortest path

• 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)

• process_obtuse (bool) – use method of Qian et al (2007) to improve accuracy for triangles with obtuse angle (default is True)

• n_secondary (int) – number of secondary nodes (SPM) (default is 5)

• n_tertiary (int) – number of tertiary nodes (DSPM) (default is 2)

• 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 2)

• tt_from_rp (bool) – compute traveltimes using raypaths (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.

static builder(filename, n_threads, cell_slowness, method, eps, maxit, process_obtuse, n_secondary, n_tertiary, radius_factor_tertiary, tt_from_rp)

Build instance of Mesh2d from VTK file

Parameters:

• filename (str) – Name of file holding a vtkUnstructuredGrid. The grid must have point or cell attribute named either ‘Slowness’, ‘slowness’, ‘Velocity’, ‘velocity’, or ‘P-wave velocity’. All cells must be of type vtkTriangle

• Constructor (Other parameters are defined in)

Returns:

mesh – mesh instance

Return type:

Mesh2d

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 (nnodes,)

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

n_threads

number of threads for raytracing

Type:

int

nparams

total number of parameters for mesh

Type:

int

raytrace(source, rcv, slowness=None, thread_no=None, aggregate_src=False, compute_L=False, return_rays=False) → tt, rays

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, (None by default)) – slowness at grid nodes or cells (depending on cell_slowness) 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 the matrix of partial derivatives of travel time w/r to the medium parameters. For a mesh with slowness defined at the cells, L holds one block of ncells columns per parameter, in the order the setters take them: slowness for an isotropic medium; slowness and xi for an elliptical one; Vs0 and gamma for vti_sh; slowness, xi and the tilt angle for tilted_elliptical; Vs0, gamma and the tilt angle for tti_sh; slowness, s2 and s4 for weakly_anelliptical; Vp0, Vs0, epsilon and delta for vti_psv; and those four with the tilt angle for tti_psv.

• 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)

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 P-wave velocity to mesh (VTI or TTI medium)

Parameters:

v (np ndarray with shape (nparams, ))

set_Vs0(v)

Assign vertical S-wave velocity to mesh (VTI or TTI medium)

Parameters:

v (np ndarray with shape (nparams, ))

set_delta(d)

Assign Thomsen’s parameter delta to mesh

Parameters:

d (np ndarray with shape (nparams, ))

set_epsilon(e)

Assign Thomsen’s parameter epsilon to mesh

Parameters:

e (np ndarray with shape (nparams, ))

set_gamma(g)

Assign Thomsen’s parameter gamma to mesh (SH wave)

Parameters:

g (np ndarray with shape (nparams, ))

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(s2)

Assign energy-velocity parameter s2 to grid

Parameters:

s2 (np ndarray with shape (nparams, ))

set_s4(s4)

Assign energy-velocity parameter s4 to grid

Parameters:

s4 (np ndarray with shape (nparams, ))

set_slowness(slowness)

Assign slowness to grid

Parameters:

slowness (np ndarray with shape (nparams, ))

set_tilt_angle(theta)

Assign tilted elliptical anisotropy angle to grid

Parameters:

theta (np ndarray with shape (nparams, ))

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 (nparams, ))

set_xi(xi)

Assign elliptical anisotropy ratio to grid

Parameters:

xi (np ndarray with shape (nparams, ))

to_vtk(fields, filename)

Save mesh 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 mesh, or a list of raypath coordinates.

• filename (str) – Name of file without extension for saving (extension vtu 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)

class ttcrpy.tmesh.Mesh3d

class to perform raytracing with tetrahedral meshes

Variables:

• nparams (int) – total number of parameters for grid

• n_threads (int) – number of threads for raytracing

• Constructor

• Mesh3d (Mesh3d(nodes, tetra, n_threads, cell_slowness, method, gradient_method, tt_from_rp, process_vel, eps, maxit, min_dist, n_secondary, n_tertiary, radius_factor_tertiary, translate_grid=False) ->) –

  param nodes:

  node coordinates

  type nodes:

  np.ndarray with shape (nnodes, 3)

  param tetra:

  indices of nodes forming the tetrahedra

  type tetra:

  np.ndarray of int with shape (ntetra, 4)

  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

  type method:

  string

  param aniso:

  type of anisotropy (SPM method and cell_slowness only)

  • ’iso’ : isotropic medium

  • ’elliptical’ : ellipsoidal anisotropy, axes aligned with the global axes; 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

  Every setter takes one value per tetrahedron. The tilted models of Mesh2d 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 gradient_method:

  method to compute traveltime gradient (default is 1)

  • 0 : least-squares first-order

  • 1 : least-squares second-order

  • 2 : Averaging-Based method

  type gradient_method:

  int

  param tt_from_rp:

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

  type tt_from_rp:

  bool

  param process_vel:

  process velocity instead of slowness at nodes when interpolating and computing matrix of partial derivative of traveltime w/r to model parameters (interpolation: for cell_slowness == False or FSM) (defauls is False)

  type process_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 min_dist:

  tolerance for backward raytracing (default is 1e-5)

  type min_dist:

  double

  param n_secondary:

  number of secondary nodes (SPM & 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

static builder(filename, n_threads, cell_slowness, method, gradient_method, tt_from_rp, process_vel, eps, maxit, min_dist, n_secondary, n_tertiary, radius_factor_tertiary, translate_grid=0)

Build instance of Mesh3d from VTK file

Parameters:

• filename (str) – Name of file holding a vtkUnstructuredGrid. The grid must have point or cell attribute named either ‘Slowness’, ‘slowness’, ‘Velocity’, ‘velocity’, or ‘P-wave velocity’. All cells must be of type vtkTetra

• Constructor (Other parameters are defined in)

Returns:

mesh – mesh instance

Return type:

Mesh3d

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)

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 (nparams,) (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(order=2, taylor_order=2, weighting=True, squared=True, s0inside=False, additional_points=0)

Compute smoothing matrices (spatial derivative)

Parameters:

• order (int) – order of derivative (1 or 2, 2 by default)

• taylor_order (int) – order of taylors series expansion (1 or 2, 2 by default)

• weighting (bool) – apply inverse distance weighting (True by default)

• squared (bool) – Second derivative evaluated by taking the square of first derivative. Applied only if order == 2 (True by default)

• s0inside (bool) – (experimental) ignore slowness value at local node (value is a filtered estimate) (False by default)

• additional_points (int) – use additional points to compute derivatives (minimum sometimes yield noisy results when rays are close to domain limits) (0 by default)

Returns:

Kx, Ky, Kz – matrices for derivatives along x, y, & z

Return type:

tuple of csr_array

data_kernel_straight_rays(Tx, Rx) → L

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

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

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 (nnodes,)

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 (nparams, ) (optional)) – slowness at grid nodes or cells (depending on cell_slowness)

Returns:

s0 – slowness at source points

Return type:

np.ndarray

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 mesh

Type:

int

raytrace(source, rcv, slowness=None, thread_no=None, aggregate_src=False, compute_L=False, return_rays=False) → tt, rays, 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, (None by default)) – slowness at grid nodes or cells (depending on cell_slowness) 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 (or velocity if process_vel == True in constructor)

• 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 (list of csr_array or scipy csr_array) – Matrix of partial derivative of travel time w/r to slowness. if input argument source has 5 columns or if slowness is defined at nodes, 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 mesh

Parameters:

v (np ndarray with shape (nparams, )) – one value per tetrahedron

set_Vs0(v)

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

Parameters:

v (np ndarray with shape (nparams, )) – one value per tetrahedron

set_chi(chi)

Assign ellipsoidal anisotropy ratio \(\chi = s_x/s_z\) to mesh

Parameters:

chi (np ndarray with shape (nparams, )) – one value per tetrahedron

set_delta(delta)

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

Parameters:

delta (np ndarray with shape (nparams, )) – one value per tetrahedron

set_epsilon(epsilon)

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

Parameters:

epsilon (np ndarray with shape (nparams, )) – one value per tetrahedron

set_gamma(gamma)

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

Parameters:

gamma (np ndarray with shape (nparams, )) – one value per tetrahedron

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 ellipsoidal anisotropy ratio \(\psi = s_y/s_z\) to mesh

Parameters:

psi (np ndarray with shape (nparams, )) – one value per tetrahedron

set_s2(s2)

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

Parameters:

s2 (np ndarray with shape (nparams, )) – one value per tetrahedron

set_s4(s4)

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

Parameters:

s4 (np ndarray with shape (nparams, )) – one value per tetrahedron

set_slowness(slowness)

Assign slowness to grid

Parameters:

slowness (np ndarray with shape (nparams, ))

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 (nparams, ))

to_vtk(fields, filename)

Save mesh 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 mesh, or a list of raypath coordinates.

• filename (str) – Name of file without extension for saving (extension vtu 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)

ttcrpy.tmesh.set_verbose(v)

Set verbosity level for C++ code

Parameters:

v (int) – verbosity level