Helper functions

JUDI provides numerous helper and utility functions need for seismic modeling and inversion.

Helper functions

Ricker wavelet
Compute CFL time stepping interval
Compute number of computational time steps
Set up 3D acquisition grid
Data interpolation
Generate and sample from frequency distribution
Read data from out of core container
Restrict model to acquisition
Additional miscellanous utilities

Ricker wavelet

Create a 1D Ricker wavelet:

Missing docstring.

Missing docstring for ricker_wavelet(tmax, dt, f0; t0=nothing). Check Documenter's build log for details.

Compute CFL time stepping interval

Calculate the time stepping interval based on the CFL condition

calculate_dt

Compute number of computational time steps

Estimate the number of computational time steps. Required for calculating the dimensions of the matrix-free linear modeling operators:

JUDI.get_computational_nt — Function

get_computational_nt(srcGeometry, recGeoemtry, model; dt=nothing)

Estimate the number of computational time steps. Required for calculating the dimensions
of the matrix-free linear modeling operators. srcGeometry and recGeometry are source
and receiver geometries of type Geometry and model is the model structure of type
Model.

get_computational_nt(Geoemtry, model; dt=nothing)

Estimate the number of computational time steps. Required for calculating the dimensions
of the matrix-free linear modeling operators. srcGeometry and recGeometry are source
and receiver geometries of type Geometry and model is the model structure of type
Model.

Set up 3D acquisition grid

Helper function to create a regular acquisition grid for a 3D survey.

setup_3D_grid

Data interpolation

Time interpolation for source/receiver data using splines. For modeling, the data is interpolated automatically onto the computational time axis, so generally, these functions are not needed for users.

JUDI.time_resample — Function

time_resample(data, geometry_in, dt_new)

Resample the input data with sinc interpolation from the current time sampling (geometrtyin) to the new time sampling `dtnew`.

Parameters

data: Data to be reampled. If data is a matrix, resamples each column.
geometry_in: Geometry on which data is defined.
dt_new: New time sampling rate to interpolate onto.

time_resample(data, dt_in, dt_new)

Resample the input data with sinc interpolation from the current time sampling dtin to the new time sampling `dtnew`.

Parameters

data: Data to be reampled. If data is a matrix, resamples each column.
dt_in: Time sampling of input
dt_new: New time sampling rate to interpolate onto.

time_resample(data, dt_in, geometry_in)

Resample the input data with sinc interpolation from the current time sampling (dtin) to the new time sampling `geometryout`.

Parameters

data: Data to be reampled. If data is a matrix, resamples each column.
geometry_out: Geometry on which data is to be interpolated.
dt_in: Time sampling rate of the data.

Generate and sample from frequency distribution

Create a probability distribution with the shape of the source spectrum from which we can draw random frequencies.

JUDI.generate_distribution — Function

generate_distribution(x; src_no=1)

Generates a probability distribution for the discrete input judiVector x.

Parameters

x: judiVector. Usualy a source with a single trace per source position.
src_no: Index of the source to select out of x

JUDI.select_frequencies — Function

select_frequencies(q_dist; fmin=0f0, fmax=Inf, nf=1)

Selects nf frequencies based on the source distribution q_dist computed with generate_distribution.

Parameters

q_dist: Distribution to sample from.
f_min: Minimum acceptable frequency to sample (defaults to 0).
f_max: Maximum acceptable frequency to sample (defaults to Inf).
fd: Number of frequnecies to sample (defaults to 1).

We can draw random samples from dist by passing it values between 0 and 1:

# Draw a single random frequency
f = dist(rand(1))

# Draw 10 random frequencies
f = dist(rand(10))

Alternatively, we can use the function:

f = select_frequencies(dist; fmin=0f0, fmax=Inf, nf=1)

to draw nf number of frequencies for a given distribution dist in the frequency range of fmin to fmax (both in kHz).

Read data from out of core container

In the case where a judiVector is out of core (points to a segy file) it is possible to convert it or part of it into an in core judiVecor with the get_data function.

d_ic = get_data(d_ooc, inds)

where inds is either a single index, a list of index or a range of index.

Restrict model to acquisition

In practice, and in particular for marine data, the aperture of a single shot is much smaller than the full model size. We provide a function (automatically used when the option limit_m is set in Options) that limits the model to the acquisition area.

JUDI.limit_model_to_receiver_area — Function

limit_model_to_receiver_area(srcGeometry, recGeometry, model, buffer; pert=nothing)

Crops the model to the area of the source an receiver with an extra buffer. This reduces the size of the problem when the model si large and the source and receiver located in a small part of the domain.

In the cartoon below, the full model will be cropped to the center area containg the source (o) receivers (x) and buffer area (*)

o Source position
x receiver positions
- Extra buffer (grid spacing in that simple case)

| . . . . . . . . . . . . . . . . . . . . . |

| . . . . . . . . . . . . . . . . . . . . . |

| . . . . * * * * * * * * * * * * . . . . . |

| . . . . * x x x x x x x x x x * . . . . . |

| . . . . * x x x x x x x x x x * . . . . . |

| . . . . * x x x x x x x x x x * . . . . . |

| . . . . * x x x x x o x x x x * . . . . . |

| . . . . * x x x x x x x x x x * . . . . . |

| . . . . * x x x x x x x x x x * . . . . . |

| . . . . * x x x x x x x x x x * . . . . . |

| . . . . * * * * * * * * * * * * . . . . . |

| . . . . . . . . . . . . . . . . . . . . . |

| . . . . . . . . . . . . . . . . . . . . . |

Parameters

srcGeometry: Geometry of the source.
recGeometry: Geometry of the receivers.
model: Model to be croped.
buffer: Size of the buffer on each side.
pert: Model perturbation (optional) to be cropped as well.

We also provide it's complement that removes receivers outside of the computational domain if the dataset contains locations that are not wanted

JUDI.remove_out_of_bounds_receivers — Function

remove_out_of_bounds_receivers(recGeometry, model)

Removes receivers that are positionned outside the computational domain defined by the model.

Parameters

recGeometry: Geometry of receivers in which out of bounds will be removed.
model: Model defining the computational domain.

remove_out_of_bounds_receivers(recGeometry, recData, model)

Removes receivers that are positionned outside the computational domain defined by the model.

Parameters

recGeometry: Geometry of receivers in which out of bounds will be removed.
recData: Shot record for that geometry in which traces will be removed.
model: Model defining the computational domain.

Additional miscellanous utilities

JUDI.devito_model — Function

devito_model(model, options;dm=nothing)

Creates a python side model strucutre for devito.

Parameters

model: JUDI Model structure.
options: JUDI Options structure.
dm: Squared slowness perturbation (optional), Array or PhysicalParameter.

JUDI.setup_grid — Function

setup_grid(geometry, n)

Sets up the coordinate arrays for Devito.

Parameters:

geometry: Geometry containing the coordinates
n: Domain size

JUDI.pad_sizes — Function

pad_sizes(model, options; so=nothing)

Computes ABC padding sizes according to the model's numbr of abc points and spatial order

Parameters

model: JUDI or Python side Model.
options: JUDI Options structure.
so: Space order (optional) defaults to options.space_order.

JUDI.pad_array — Function

pad_array(m, nb; mode=:border)

Pads to the input array with either copying the edge value (:border) or zeros (:zeros)

Parameters

m: Array to be padded.
nb: Size of padding. Array of tuple with one (nbleft, nbright) tuple per dimension.
mode: Padding mode (optional), defaults to :border.

JUDI.remove_padding — Function

remove_padding(m, nb; true_adjoint=False)

Removes the padding from array m. This is the adjoint of pad_array.

Parameters

m: Array to remvove padding from.
nb: Size of padding. Array of tuple with one (nbleft, nbright) tuple per dimension.
true_adjoint: Unpadding mode, defaults to False. Will sum the padding to the edge point with true_adjoint=true

and should only be used this way for adjoint testing purpose.

JUDI.convertToCell — Function

convertToCell(x)

Convert an array x to a cell array (Array{Any,1}) with length(x) entries,
where the i-th cell contains the i-th entry of x.

Parameters

x: Array to be converted into and array of array

JUDI.process_input_data — Function

process_input_data(input, geometry, nsrc)

Preprocesses input Array into an Array of Array for modeling

Parameters:

input: Input to preprocess.
geometry: Geometry containing physical parameters.
nsrc: Number of sources

process_input_data(input, model, nsrc)

Preprocesses input Array into an Array of Array for modeling

Parameters:

input: Input to preprocess.
model: Model containing physical parameters.
nsrc: Number of sources

Base.reshape — Function

reshape(A, dims...) -> AbstractArray
reshape(A, dims) -> AbstractArray

Return an array with the same data as A, but with different dimension sizes or number of dimensions. The two arrays share the same underlying data, so that the result is mutable if and only if A is mutable, and setting elements of one alters the values of the other.

The new dimensions may be specified either as a list of arguments or as a shape tuple. At most one dimension may be specified with a :, in which case its length is computed such that its product with all the specified dimensions is equal to the length of the original array A. The total number of elements must not change.

Examples

julia> A = Vector(1:16)
16-element Vector{Int64}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16

julia> reshape(A, (4, 4))
4×4 Matrix{Int64}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> reshape(A, 2, :)
2×8 Matrix{Int64}:
 1  3  5  7   9  11  13  15
 2  4  6  8  10  12  14  16

julia> reshape(1:6, 2, 3)
2×3 reshape(::UnitRange{Int64}, 2, 3) with eltype Int64:
 1  3  5
 2  4  6

JUDI.transducer — Function

transducer(q, d, r, theta)

Create the JUDI soure for a circular transducer Theta=0 points downward:

. . . . - - - . . . . . .

. . . . + + + . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

Theta=pi/2 points right:

. . . . - + . . . . . . .

. . . . - + . . . . . . .

. . . . - + . . . . . . .

. . . . . . . . . . . . .

2D only, to extend to 3D

JUDI.as_vec — Function

as_vec(x, ::Val{Bool})

Vectorizes output when return_array is set to true.