Edit on GitHub

Reference

JetPackWaveFD.JopLnProp2DAcoIsoDenQ_DEO2_FDTD — Method
JopNlProp2DAcoIsoDenQ_DEO2_FDTD(; kwargs...)
JopLnProp2DAcoIsoDenQ_DEO2_FDTD(; v, kwargs...)

Create a Jets nonlinear or linearized operator for 2D visco-acoustic, isotropic, variable density modeling.

Model Parameters

This propagator operates with two model parameters, as shown in the table below.

ParameterDescription
vP wave velocity
bbuoyancy (reciprocal density)

Velocity v and buoyancy b can be active parameters and inverted for using the Jacobian machinery, or passive parameters that are constant. If a parameter is active there will be a component of model assigned to it.

The model is 3D with the slowest dimension for model parameter. The earth model property maps to the slow- dimension array index via the state(F).active_modelset dictionary.

Examples

Model and acquisition geometry setup

  1. load modules Jets, WaveFD, and JetPackWaveFD
  2. set up the model discretization, coordinate size and spacing
  3. set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
  4. create constant buoyancy array
using Jets, WaveFD, JetPackWaveFD
nz,nx = 200, 100                      # spatial discretization size
dz,dx = 20.0, 20.0                    # spatial discretization sampling
ntrec = 1101                          # number of temporal samples in recorded data
dtrec = 0.0100                        # sample rate for recorded data
dtmod = 0.0025                        # sample rate for modeled data
wavelet = WaveletCausalRicker(f=5.0)  # type of wavelet to use (source signature) 
sz = dz                               # source Z location 
sx = dx*(nx/2)                        # source X location 
rx = dx*[0:0.5:nx-1;];                # array of receiver X locations
rz = 2*dz*ones(length(0:0.5:nx-1));   # array of receiver Z locations
b = ones(Float32, nz, nx);            # buoyancy model (reciprocal density)

Construct and apply the nonlinear operator

  1. create the nonlinear operator F
  2. create the constant velocity model m₀
  3. perform nonlinear forward modeling with constant velocity model m₀ and return the resulting modeled data in d
F = JopNlProp2DAcoIsoDenQ_DEO2_FDTD(; b = b, ntrec = ntrec, 
    dz = dz, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = 1500 .* ones(domain(F));
d = F*m₀;              # forward nonlinear op

Construct and apply the linearized Jacobian operator (method 1)

  1. create the nonlinear operator F directly by construction of a jet.
  2. create the constant velocity model m₀
  3. create the Jacobian operator J by linearization of F at point m₀
  4. create a random model perturbation vector δm.
  5. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  6. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm.

Note that the Jacobian operators J and J' require the serialized nonlinear forward wavefield, and this is generated automatically whenever required. If you watch the logging to standard out from this example, you will first see the finite difference evolution for the nonlinear forward, followed by the linearized forward, and finally the linearized adjoint.

F = JopNlProp2DAcoIsoDenQ_DEO2_FDTD(; b = b, ntrec = ntrec, 
    dz = dz, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = 1500 .* ones(domain(F));
J = jacobian(F, m₀)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Construct and apply the linearized Jacobian operator (method 2)

  1. create the constant velocity model m₀
  2. create the Jacobian operator J at point m₀ directly by construction of a jet.
  3. create a random model perturbation vector δm.
  4. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  5. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm. 
m₀ = 1500 * ones(Float32, nz, nx);
J = JopLnProp2DAcoIsoDenQ_DEO2_FDTD(; v = m₀, b = b, ntrec = ntrec, 
    dz = dz, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sx = sx, rz = rz, rx = rx)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Required Parameters

  • v the point at which the jet is linearized. Note this argument is required in the constuctor for JopLnProp2DAcoIsoDenQ_DEO2_FDTD but not JopNlProp2DAcoIsoDenQ_DEO2_FDTD. This constuctor is shown in the last example above.
  • dtmod The sample rate for the modeled data. You can establish a lower bound for the modeling sample rate with the following expression: dt = 0.75 * 0.38 * max(dx,dz) / maximum(v). Note that the usual Courant–Friedrichs–Lewy condition (CFL condition) for stability in finite difference modeling is modified heuristically to include the impact of the visco-acoustic implementation of this operator, requiring a 25% smaller smaple rate.
  • dtrec The number of time samples in the recorded data. Note requirement that dtrec be an even multiple of dtmod.
  • ntrec The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined from ntrec, dtrec, and dtmod.

Optional Parameters

Defaults for arguments are shown inside square brackets.

  • b [ones(Float32,0,0)] the buoyancy (reciprocal density) array.
  • srcfieldfile [joinpath(tempdir(), "field-c3628917-2e59-43dc-9fb7-c58a2f8e1d4d.bin")] the full path to a scratch file used for the serializationof the compressed nonlinear source wavefield.
  • comptype [nothing] the type of compression to use for the serialization of the nonlinear source wavefield. The type of compression must be one of:
    • nothing - no compression.
    • Float32 - if wavefield is Float64, then do a simple conversion to Float32.
    • UInt32 - compression using CvxCompress (windowing + 2D wavelet transform + thresholding + quantization + run-length-encoding).
  • compscale [1e-2] determines the thresholding for the compression of the nonlinear source wavefield prior to serialization. Larger values mean more aggressive compression. You can likely increase compscale to 1.0 before you start to notice differences in output from Jacobian operations.
  • nz_subcube, nx_subcube [32] The Z and X sizes of windows used for compression of the nonlinear source wavefield with CvxCompress. Note the requirement [8 <= n*_subcube <=256].
  • isinterior [false] boolean flag that indicates how the nonlinear source wavefield is saved. For large models, operation will be faster with isinterior = true, but the linearization correctness test may fail.
    • true the entire model including absorbing boundaries is serialized and deserialized
    • false the interior part of the model excluding absorbing boundaries is serialized and deserialized
  • sz [0.0] Array of source Z coordinate. Note that if multiple sources are provided, they will will be injected simultaneously during finite difference evolution.
  • sx [0.0] Array of source X coordinate.
  • st [0.0] Array of source delay times.
  • interpmethod [:hicks] Type of physical interpolation for sources and receivers. For locations that are not on the physical grid coordinates, interpolation is used either to inject or extract information. interpmethod must be one of:
    • :hicks Hicks 2D sinc interpolation (up to 8x8 nonzero points per location)
    • :linear bilinear interpolation (up to 2x2 nonzero points per location)
  • rz [[0.0]] 2D array of receiver Z coordinates
  • rx [[0.0]] 2D array of receiver Z coordinates
  • z0 [0.0] Origin of physical coordinates for the Z dimension.
  • x0 [0.0] Origin of physical coordinates for the X dimension.
  • dz [10.0] Spacing of physical coordinates in the Z dimension.
  • dx [10.0] Spacing of physical coordinates in the X dimension.
  • freqQ [5.0] The center frequency for the Maxwell body approximation to dissipation only attenuation. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qMin [0.1] The minimum value for Qp at the boundary of the model used in our Maxwell body approximation to dissipation only attenuation. This is not a physically meaningful value for Qp, as we use the attenuation to implement absorbing boundary conditions and eliminate outgoing waves on the boundaries of the computational domain. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qInterior [100.0] the value for Qp in the interior of the model used in our Maxwell body approximation to dissipation only attenuation. This is the value for Qp away from the absorbing boundaries and is a physically meaningful value. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • padz, padx [0.0], [0.0] - apply extra padding to the survey determined aperture in Ginsu. Please see Ginsu for more information.
  • nbz_cache, nbx_cache [512], [8] The size of cache blocks in the Z and X dimensions. In general the cache block in the Z (fast) dimension should be ≥ the entire size of that dimension, and the cache block size in the slower dimension is generally small in order to allow the entire block to fit in cache.
  • nsponge [50] The number of grid cells to use for the absorbing boundary. For high fidelity modeling this should be > 60 grid points, but can be significantly smaller for some use cases like low frequency full waveform inversion.
  • wavelet [WaveletCausalRicker(f=5.0)] The source wavelet, can be specified as either a Wavelet type or an array.
  • freesurface [false] Determines if a free surface (true) or absorbing (false) top boundary condition is applied.
  • imgcondition ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint only exists for "standard" imaging condition currently.
  • RTM_weight [0.5] determines the balance of short wavelengths and long wavelengths in the imaging condition. A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition.
  • nthreads [Sys.CPU_THREADS] The number of threads to use for OpenMP parallelization of the modeling.
  • reportinterval [500] The interval at which information about the propagtion is logged.

See also: Ginsu, WaveletSine, WaveletRicker, WaveletMinPhaseRicker, WaveletDerivRicker, WaveletCausalRicker, WaveletOrmsby, WaveletMinPhaseOrmsby

source
JetPackWaveFD.JopLnProp2DAcoTTIDenQ_DEO2_FDTD — Method
JopNlProp2DAcoTTIDenQ_DEO2_FDTD(; kwargs...)
JopLnProp2DAcoTTIDenQ_DEO2_FDTD(; m₀, kwargs...)

Create a Jets nonlinear or linearized operator for 2D pseudo- visco-acoustic, tilted transverse isotropy, variable density modeling.

Model Parameters

This propagator operates with five model parameters, as shown in the table below.

ParameterDescription
vP wave velocity
ϵModified Thomsen's weak anisotropy parameter
ηModified Alkhalifah's weak anisotropy parameter η
bbuoyancy (reciprocal density)
θsymmetry axis tilt angle from vertical (radians)

Pseudo-acoustic approximation

With the pseudo-acoustic approximation we employ, the transformation to self-adjoint form requires a parameter f representing the average ratio of shear-wave to p-wave velocity, and a modfication to Alkhalifah's weak anisotropy parameter η. We show formulas for these two parameters below. The parameter f is specified as a scalar for the entire model, with default value 0.85, implying a shear velocity around 38% of the p wave velocity.

  • η = sqrt[2(ϵ-δ)/(f+2ϵ)]
  • f = 1 - Vₛ² / Vₚ²

For more information about the pseudo-acoustic approximation we employ please see https://library.seg.org/doi/10.1190/segam2016-13878451.1.

Active and passive model parameters

There are two modes of operation that define different sets of active and passive parameters, velocity only and velocity and anisotropy.

An active parameter can be inverted for using the Jacobian linearized operator machinery, and a passive parameter is constant. The two modes are:

ModeActive ParametersPassive Parameters
velocity only[v][ϵ, η, b, θ ]
velocity and anisotropy[v, ϵ, η][b, θ ]

To make a parameter passive, you pass the array for that parameter in the constructor for the operator, implying that it is in state and does not change with the action of the operator. Parameters that are not passed to the constructor are assumed to be active, and will be part of the model that the operator acts on, stored as a 3D array with the following indices: [Z coord][X coord][Active Parameter]

Examples

Model and acquisition geometry setup

  1. load modules Jets, WaveFD, and JetPackWaveFD
  2. set up the model discretization, coordinate size and spacing
  3. set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
  4. create constant models for v, ϵ, η, b, θ
using Jets, WaveFD, JetPackWaveFD
nz,nx = 200, 100                      # spatial discretization size
dz,dx = 20.0, 20.0                    # spatial discretization sampling
ntrec = 1101                          # number of temporal samples in recorded data
dtrec = 0.0100                        # sample rate for recorded data
dtmod = 0.0025                        # sample rate for modeled data
wavelet = WaveletCausalRicker(f=5.0)  # type of wavelet to use (source signature) 
sz = dz                               # source Z location 
sx = dx*(nx/2)                        # source X location 
rx = dx*[0:0.5:nx-1;];                # array of receiver X locations
rz = 2*dz*ones(length(0:0.5:nx-1));   # array of receiver Z locations
b = ones(Float32, nz, nx);            # constant buoyancy
ϵ = 0.1 * ones(Float32, nz, nx);      # constant epsilon
η = 0.2 * ones(Float32, nz, nx);      # constant eta
θ = (π/8) * ones(Float32, nz, nx);    # constant tilt angle

Construct and apply the nonlinear operator (v is active; ϵ, η, b, θ are passive)

  1. create the nonlinear operator F
  2. create the model vector and set parameter v
  3. perform nonlinear forward modeling with model m₀ and return the resulting modeled data in d
F = JopNlProp2DAcoTTIDenQ_DEO2_FDTD(; b = b, ϵ = ϵ, η = η, θ = θ, ntrec = ntrec, dtrec = dtrec, 
    dtmod = dtmod, dz = dz, dx = dx, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,1] .= 1500;
d = F*m₀;              # forward nonlinear op

Construct and apply the nonlinear operator (v, ϵ, η are active; b, θ are passive)

  1. create the nonlinear operator F
  2. create the model vector and set parameters v, ϵ, η
  3. perform nonlinear forward modeling with model m₀ and return the resulting modeled data in d
F = JopNlProp2DAcoTTIDenQ_DEO2_FDTD(; b = b, θ = θ, ntrec = ntrec, dtrec = dtrec, dtmod = dtmod, 
    dz = dz, dx = dx, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,1] .= 1500;
m₀[:,:,2] .= 0.1;
m₀[:,:,3] .= 0.2;
d = F*m₀;              # forward nonlinear op

Construct and apply the linearized Jacobian operator (method 1)

For this example we assume in the model that v, ϵ, η are active parameters, and b, θ are passive.

  1. create the nonlinear operator F directly by construction of a jet.
  2. create the model vector and set parameters v, ϵ, η
  3. create the Jacobian operator J by linearization of F at point m₀
  4. create a random model perturbation vector δm.
  5. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  6. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm.

Note that the Jacobian operators J and J' require the serialized nonlinear forward wavefield, and this is generated automatically whenever required. If you watch the logging to standard out from this example, you will first see the finite difference evolution for the nonlinear forward, followed by the linearized forward, and finally the linearized adjoint.

F = JopNlProp2DAcoTTIDenQ_DEO2_FDTD(; b = b, θ = θ, ntrec = ntrec, dz = dz, dx = dx, 
    dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,1] .= 1500;
m₀[:,:,2] .= 0.1;
m₀[:,:,3] .= 0.2;
J = jacobian(F, m₀)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Construct and apply the linearized Jacobian operator (method 2)

For this example we assume in the model that v, ϵ, η are active parameters, and b, θ are passive.

  1. create the model vector and set parameters v, ϵ, η
  2. create the Jacobian operator J at point m₀ directly by construction of a jet.
  3. create a random model perturbation vector δm.
  4. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  5. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm. 
m₀ = zeros(Float32, nz, nx, 3);
m₀[:,:,1] .= 1500;
m₀[:,:,2] .= 0.1;
m₀[:,:,3] .= 0.2;
J = JopLnProp2DAcoTTIDenQ_DEO2_FDTD(; m₀ = m₀, b = b, θ = θ, ntrec = ntrec, dz = dz, dx = dx, 
    dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Required Parameters

  • m₀ the point at which the jet is linearized. Note this argument is required in the constuctor for JopLnProp2DAcoTTIDenQ_DEO2_FDTD but not JopNlProp2DAcoTTIDenQ_DEO2_FDTD. This constuctor is shown in the last example above. Please note that you must consider which parameters are active and passive, per the discussion on model parameters above and examples.
  • dtmod The sample rate for the modeled data. You can establish a lower bound for the modeling sample rate with the following expression: dt = 0.75 * 0.38 * max(dx,dz) / maximum(v). Note that the usual Courant–Friedrichs–Lewy condition (CFL condition) for stability in finite difference modeling is modified heuristically to include the impact of the visco-acoustic implementation of this operator, requiring a 25% smaller smaple rate.
  • dtrec The number of time samples in the recorded data. Note requirement that dtrec be an even multiple of dtmod.
  • ntrec The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined from ntrec, dtrec, and dtmod.

Optional Parameters

Defaults for arguments are shown inside square brackets.

  • b [ones(Float32,0,0)] the buoyancy (reciprocal density) array.
  • θ [zeros(Float32,0,0)] the symmetry axies tilt angle in radians. Assumed to be zero if not specified. θ must be an array the same size as buoyancy b.
  • f [0.85] See the discussion above regarding the Pseudo-acoustic approximation used here
  • srcfieldfile [joinpath(tempdir(), "field-55e40154-73c6-434c-b029-1fe5243e0e7d.bin")] the full path to a scratch file used for the serializationof the compressed nonlinear source wavefield.
  • comptype [nothing] the type of compression to use for the serialization of the nonlinear source wavefield. The type of compression must be one of:
    • nothing - no compression.
    • Float32 - if wavefield is Float64, then do a simple conversion to Float32.
    • UInt32 - compression using CvxCompress (windowing + 2D wavelet transform + thresholding + quantization + run-length-encoding).
  • compscale [1e-2] determines the thresholding for the compression of the nonlinear source wavefield prior to serialization. Larger values mean more aggressive compression. You can likely increase compscale to 1.0 before you start to notice differences in output from Jacobian operations.
  • nz_subcube, nx_subcube [32] The Z and X sizes of windows used for compression of the nonlinear source wavefield with CvxCompress. Note the requirement [8 <= n*_subcube <=256].
  • isinterior [false] boolean flag that indicates how the nonlinear source wavefield is saved. For large models, operation will be faster with isinterior = true, but the linearization correctness test may fail.
    • true the entire model including absorbing boundaries is serialized and deserialized
    • false the interior part of the model excluding absorbing boundaries is serialized and deserialized
  • sz [0.0] Array of source Z coordinate. Note that if multiple sources are provided, they will will be injected simultaneously during finite difference evolution.
  • sx [0.0] Array of source X coordinate.
  • st [0.0] Array of source delay times.
  • interpmethod [:hicks] Type of physical interpolation for sources and receivers. For locations that are not on the physical grid coordinates, interpolation is used either to inject or extract information. interpmethod must be one of:
    • :hicks Hicks 2D sinc interpolation (up to 8x8 nonzero points per location)
    • :linear bilinear interpolation (up to 2x2 nonzero points per location)
  • rz [[0.0]] 2D array of receiver Z coordinates
  • rx [[0.0]] 2D array of receiver Z coordinates
  • z0 [0.0] Origin of physical coordinates for the Z dimension.
  • x0 [0.0] Origin of physical coordinates for the X dimension.
  • dz [10.0] Spacing of physical coordinates in the Z dimension.
  • dx [10.0] Spacing of physical coordinates in the X dimension.
  • freqQ [5.0] The center frequency for the Maxwell body approximation to dissipation only attenuation. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qMin [0.1] The minimum value for Qp at the boundary of the model used in our Maxwell body approximation to dissipation only attenuation. This is not a physically meaningful value for Qp, as we use the attenuation to implement absorbing boundary conditions and eliminate outgoing waves on the boundaries of the computational domain. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qInterior [100.0] the value for Qp in the interior of the model used in our Maxwell body approximation to dissipation only attenuation. This is the value for Qp away from the absorbing boundaries and is a physically meaningful value. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • padz, padx [0.0], [0.0] - apply extra padding to the survey determined aperture in Ginsu. Please see Ginsu for more information.
  • nbz_cache, nbx_cache [512], [8] The size of cache blocks in the Z and X dimensions. In general the cache block in the Z (fast) dimension should be ≥ the entire size of that dimension, and the cache block size in the slower dimension is generally small in order to allow the entire block to fit in cache.
  • nsponge [50] The number of grid cells to use for the absorbing boundary. For high fidelity modeling this should be > 60 grid points, but can be significantly smaller for some use cases like low frequency full waveform inversion.
  • wavelet [WaveletCausalRicker(f=5.0)] The source wavelet, can be specified as either a Wavelet type or an array.
  • freesurface [false] Determines if a free surface (true) or absorbing (false) top boundary condition is applied.
  • imgcondition ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint only exists for "standard" imaging condition currently.
  • RTM_weight [0.5] determines the balance of short wavelengths and long wavelengths in the imaging condition. A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition.
  • nthreads [Sys.CPU_THREADS] The number of threads to use for OpenMP parallelization of the modeling.
  • reportinterval [500] The interval at which information about the propagtion is logged.

See also: Ginsu, WaveletSine, WaveletRicker, WaveletMinPhaseRicker, WaveletDerivRicker, WaveletCausalRicker, WaveletOrmsby, WaveletMinPhaseOrmsby

source
JetPackWaveFD.JopLnProp2DAcoVTIDenQ_DEO2_FDTD — Method
JopNlProp2DAcoVTIDenQ_DEO2_FDTD(; kwargs...)
JopLnProp2DAcoVTIDenQ_DEO2_FDTD(; m₀, kwargs...)

Create a Jets nonlinear or linearized operator for 2D pseudo- visco-acoustic, vertical transverse isotropy, variable density modeling.

Model Parameters

This propagator operates with four model parameters, as shown in the table below.

ParameterDescription
vP wave velocity
ϵModified Thomsen's weak anisotropy parameter
ηModified Alkhalifah's weak anisotropy parameter η
bbuoyancy (reciprocal density)

Pseudo-acoustic approximation

With the pseudo-acoustic approximation we employ, the transformation to self-adjoint form requires a parameter f representing the average ratio of shear-wave to p-wave velocity, and a modfication to Alkhalifah's weak anisotropy parameter η. We show formulas for these two parameters below. The parameter f is specified as a scalar for the entire model, with default value 0.85, implying a shear velocity around 38% of the p wave velocity.

  • η = sqrt[2(ϵ-δ)/(f+2ϵ)]
  • f = 1 - Vₛ² / Vₚ²

For more information about the pseudo-acoustic approximation we employ please see https://library.seg.org/doi/10.1190/segam2016-13878451.1.

Active and passive model parameters

There are two modes of operation that define different sets of active and passive parameters, velocity only and velocity and anisotropy.

An active parameter can be inverted for using the Jacobian linearized operator machinery, and a passive parameter is constant. The two modes are:

ModeActive ParametersPassive Parameters
velocity only[v][ϵ, η, b]
velocity and anisotropy[v, ϵ, η][b]

To make a parameter passive, you pass the array for that parameter in the constructor for the operator, implying that it is in state and does not change with the action of the operator. Parameters that are not passed to the constructor are assumed to be active, and will be part of the model that the operator acts on, stored as a 3D array with the following indices: [Z coord][X coord][Active Parameter]

Examples

Model and acquisition geometry setup

  1. load modules Jets, WaveFD, and JetPackWaveFD
  2. set up the model discretization, coordinate size and spacing
  3. set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
  4. create constant models for v, ϵ, η, b
using Jets, WaveFD, JetPackWaveFD
nz,nx = 200, 100                      # spatial discretization size
dz,dx = 20.0, 20.0                    # spatial discretization sampling
ntrec = 1101                          # number of temporal samples in recorded data
dtrec = 0.0100                        # sample rate for recorded data
dtmod = 0.0025                        # sample rate for modeled data
wavelet = WaveletCausalRicker(f=5.0)  # type of wavelet to use (source signature) 
sz = dz                               # source Z location 
sx = dx*(nx/2)                        # source X location 
rx = dx*[0:0.5:nx-1;];                # array of receiver X locations
rz = 2*dz*ones(length(0:0.5:nx-1));   # array of receiver Z locations
b = ones(Float32, nz, nx);            # constant buoyancy
ϵ = 0.1 * ones(Float32, nz, nx);      # constant epsilon
η = 0.2 * ones(Float32, nz, nx);      # constant eta

Construct and apply the nonlinear operator (v is active; ϵ, η, b are passive)

  1. create the nonlinear operator F
  2. create the model vector and set parameter v
  3. perform nonlinear forward modeling with model m₀ and return the resulting modeled data in d
F = JopNlProp2DAcoVTIDenQ_DEO2_FDTD(; b = b, ϵ = ϵ, η = η, ntrec = ntrec, dtrec = dtrec, 
    dtmod = dtmod, dz = dz, dx = dx, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,1] .= 1500;
d = F*m₀;              # forward nonlinear op

Construct and apply the nonlinear operator (v, ϵ, η are active; b is passive)

  1. create the nonlinear operator F
  2. create the model vector and set parameters v, ϵ, η
  3. perform nonlinear forward modeling with model m₀ and return the resulting modeled data in d
F = JopNlProp2DAcoVTIDenQ_DEO2_FDTD(; b = b, ntrec = ntrec, dtrec = dtrec, dtmod = dtmod,
    dz = dz, dx = dx, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,1] .= 1500;
m₀[:,:,2] .= 0.1;
m₀[:,:,3] .= 0.2;
d = F*m₀;              # forward nonlinear op

Construct and apply the linearized Jacobian operator (method 1)

For this example we assume in the model that v, ϵ, η are active parameters, and b is passive.

  1. create the nonlinear operator F directly by construction of a jet.
  2. create the model vector and set parameters v, ϵ, η
  3. create the Jacobian operator J by linearization of F at point m₀
  4. create a random model perturbation vector δm.
  5. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  6. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm.

Note that the Jacobian operators J and J' require the serialized nonlinear forward wavefield, and this is generated automatically whenever required. If you watch the logging to standard out from this example, you will first see the finite difference evolution for the nonlinear forward, followed by the linearized forward, and finally the linearized adjoint.

F = JopNlProp2DAcoVTIDenQ_DEO2_FDTD(; b = b, ntrec = ntrec, dz = dz, dx = dx, 
    dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,1] .= 1500;
m₀[:,:,2] .= 0.1;
m₀[:,:,3] .= 0.2;
J = jacobian(F, m₀)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Construct and apply the linearized Jacobian operator (method 2)

For this example we assume in the model that v, ϵ, η are active parameters, and b is passive.

  1. create the model vector and set parameters v, ϵ, η
  2. create the Jacobian operator J at point m₀ directly by construction of a jet.
  3. create a random model perturbation vector δm.
  4. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  5. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm. 
m₀ = zeros(Float32, nz, nx, 3);
m₀[:,:,1] .= 1500;
m₀[:,:,2] .= 0.1;
m₀[:,:,3] .= 0.2;
J = JopLnProp2DAcoVTIDenQ_DEO2_FDTD(; m₀ = m₀, b = b, ntrec = ntrec, dz = dz, dx = dx, 
    dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Required Parameters

  • m₀ the point at which the jet is linearized. Note this argument is required in the constuctor for JopLnProp2DAcoVTIDenQ_DEO2_FDTD but not JopNlProp2DAcoVTIDenQ_DEO2_FDTD. This constuctor is shown in the last example above. Please note that you must consider which parameters are active and passive, per the discussion on model parameters above and examples.
  • dtmod The sample rate for the modeled data. You can establish a lower bound for the modeling sample rate with the following expression: dt = 0.75 * 0.38 * max(dx,dz) / maximum(v). Note that the usual Courant–Friedrichs–Lewy condition (CFL condition) for stability in finite difference modeling is modified heuristically to include the impact of the visco-acoustic implementation of this operator, requiring a 25% smaller smaple rate.
  • dtrec The number of time samples in the recorded data. Note requirement that dtrec be an even multiple of dtmod.
  • ntrec The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined from ntrec, dtrec, and dtmod.

Optional Parameters

Defaults for arguments are shown inside square brackets.

  • b [ones(Float32,0,0)] the buoyancy (reciprocal density) array.
  • f [0.85] See the discussion above regarding the Pseudo-acoustic approximation used here
  • srcfieldfile [joinpath(tempdir(), "field-24b82d11-9647-4676-a2e8-c34c19b03b99.bin")] the full path to a scratch file used for the serializationof the compressed nonlinear source wavefield.
  • comptype [nothing] the type of compression to use for the serialization of the nonlinear source wavefield. The type of compression must be one of:
    • nothing - no compression.
    • Float32 - if wavefield is Float64, then do a simple conversion to Float32.
    • UInt32 - compression using CvxCompress (windowing + 2D wavelet transform + thresholding + quantization + run-length-encoding).
  • compscale [1e-2] determines the thresholding for the compression of the nonlinear source wavefield prior to serialization. Larger values mean more aggressive compression. You can likely increase compscale to 1.0 before you start to notice differences in output from Jacobian operations.
  • nz_subcube, nx_subcube [32] The Z and X sizes of windows used for compression of the nonlinear source wavefield with CvxCompress. Note the requirement [8 <= n*_subcube <=256].
  • isinterior [false] boolean flag that indicates how the nonlinear source wavefield is saved. For large models, operation will be faster with isinterior = true, but the linearization correctness test may fail.
    • true the entire model including absorbing boundaries is serialized and deserialized
    • false the interior part of the model excluding absorbing boundaries is serialized and deserialized
  • sz [0.0] Array of source Z coordinate. Note that if multiple sources are provided, they will will be injected simultaneously during finite difference evolution.
  • sx [0.0] Array of source X coordinate.
  • st [0.0] Array of source delay times.
  • interpmethod [:hicks] Type of physical interpolation for sources and receivers. For locations that are not on the physical grid coordinates, interpolation is used either to inject or extract information. interpmethod must be one of:
    • :hicks Hicks 2D sinc interpolation (up to 8x8 nonzero points per location)
    • :linear bilinear interpolation (up to 2x2 nonzero points per location)
  • rz [[0.0]] 2D array of receiver Z coordinates
  • rx [[0.0]] 2D array of receiver Z coordinates
  • z0 [0.0] Origin of physical coordinates for the Z dimension.
  • x0 [0.0] Origin of physical coordinates for the X dimension.
  • dz [10.0] Spacing of physical coordinates in the Z dimension.
  • dx [10.0] Spacing of physical coordinates in the X dimension.
  • freqQ [5.0] The center frequency for the Maxwell body approximation to dissipation only attenuation. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qMin [0.1] The minimum value for Qp at the boundary of the model used in our Maxwell body approximation to dissipation only attenuation. This is not a physically meaningful value for Qp, as we use the attenuation to implement absorbing boundary conditions and eliminate outgoing waves on the boundaries of the computational domain. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qInterior [100.0] the value for Qp in the interior of the model used in our Maxwell body approximation to dissipation only attenuation. This is the value for Qp away from the absorbing boundaries and is a physically meaningful value. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • padz, padx [0.0], [0.0] - apply extra padding to the survey determined aperture in Ginsu. Please see Ginsu for more information.
  • nbz_cache, nbx_cache [512], [8] The size of cache blocks in the Z and X dimensions. In general the cache block in the Z (fast) dimension should be ≥ the entire size of that dimension, and the cache block size in the slower dimension is generally small in order to allow the entire block to fit in cache.
  • nsponge [50] The number of grid cells to use for the absorbing boundary. For high fidelity modeling this should be > 60 grid points, but can be significantly smaller for some use cases like low frequency full waveform inversion.
  • wavelet [WaveletCausalRicker(f=5.0)] The source wavelet, can be specified as either a Wavelet type or an array.
  • freesurface [false] Determines if a free surface (true) or absorbing (false) top boundary condition is applied.
  • imgcondition ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint only exists for "standard" imaging condition currently.
  • RTM_weight [0.5] determines the balance of short wavelengths and long wavelengths in the imaging condition. A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition.
  • nthreads [Sys.CPU_THREADS] The number of threads to use for OpenMP parallelization of the modeling.
  • reportinterval [500] The interval at which information about the propagtion is logged.

See also: Ginsu, WaveletSine, WaveletRicker, WaveletMinPhaseRicker, WaveletDerivRicker, WaveletCausalRicker, WaveletOrmsby, WaveletMinPhaseOrmsby

source
JetPackWaveFD.JopLnProp3DAcoIsoDenQ_DEO2_FDTD — Method
JopNlProp3DAcoIsoDenQ_DEO2_FDTD(; kwargs...)
JopLnProp3DAcoIsoDenQ_DEO2_FDTD(; v, kwargs...)

Create a Jets nonlinear or linearized operator for 3D visco-acoustic, isotropic, variable density modeling.

Model Parameters

This propagator operates with two model parameters, as shown in the table below.

ParameterDescription
vP wave velocity
bbuoyancy (reciprocal density)

Velocity v and buoyancy b can be active parameters and inverted for using the Jacobian machinery, or passive parameters that are constant. If a parameter is active there will be a component of model assigned to it.

The model is 4D with the slowest dimension for model parameter. The earth model property maps to the slow- fimension array index via the state(F).active_modelset dictionary.

Examples

Model and acquisition geometry setup

  1. load modules Jets, WaveFD, and JetPackWaveFD
  2. set up the model discretization, coordinate size and spacing
  3. set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
  4. create constant buoyancy array
using Jets, WaveFD, JetPackWaveFD
nz,ny,nx = 100, 80, 60                        # spatial discretization size
dz,dy,dx = 20.0, 20.0, 20.0                   # spatial discretization sampling
ntrec = 1101                                  # number of temporal samples in recorded data
dtrec = 0.0100                                # sample rate for recorded data
dtmod = 0.0025                                # sample rate for modeled data
wavelet = WaveletCausalRicker(f=5.0)          # type of wavelet to use (source signature) 
sz = dz                                       # source Z location 
sy = dy*(ny/2)                                # source Y location 
sx = dx*(nx/2)                                # source X location 
rz = [dz for iy = 1:ny, ix = 1:nx][:];        # Array of receiver Y locations 
ry = [(iy-1)*dy for iy = 1:ny, ix = 1:nx][:]; # Array of receiver Y locations
rx = [(ix-1)*dx for iy = 1:ny, ix=1:nx][:];   # Array of receiver X locations
b = ones(Float32, nz, ny, nx);                # buoyancy model (reciprocal density)

Construct and apply the nonlinear operator

  1. create the nonlinear operator F
  2. create the constant velocity model m₀
  3. perform nonlinear forward modeling with constant velocity model m₀ and return the resulting modeled data in d
F = JopNlProp3DAcoIsoDenQ_DEO2_FDTD(; b = b, isinterior=true, nsponge = 10, ntrec = ntrec, 
    dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = 1500 .* ones(domain(F));
d = F*m₀;              # forward nonlinear op

Construct and apply the linearized Jacobian operator (method 1)

  1. create the nonlinear operator F directly by construction of a jet.
  2. create the constant velocity model m₀
  3. create the Jacobian operator J by linearization of F at point m₀
  4. create a random model perturbation vector δm.
  5. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  6. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm.

Note that the Jacobian operators J and J' require the serialized nonlinear forward wavefield, and this is generated automatically whenever required. If you watch the logging to standard out from this example, you will first see the finite difference evolution for the nonlinear forward, followed by the linearized forward, and finally the linearized adjoint.

F = JopNlProp3DAcoIsoDenQ_DEO2_FDTD(; b = b, isinterior=true, nsponge = 10, ntrec = ntrec, 
    dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = 1500 .* ones(domain(F));
J = jacobian(F, m₀)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Construct and apply the linearized Jacobian operator (method 2)

  1. create the constant velocity model m₀
  2. create the Jacobian operator J at point m₀ directly by construction of a jet.
  3. create a random model perturbation vector δm.
  4. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  5. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm. 
m₀ = 1500 .* ones(Float32, nz, ny, nx);       # velocity model
J = JopLnProp3DAcoIsoDenQ_DEO2_FDTD(; v = m₀, b = b, isinterior=true, nsponge = 10, ntrec = ntrec, 
    dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sy = sy, sx = sx, ry = ry, rz = rz, rx = rx)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Required Parameters

  • v the point at which the jet is linearized. Note this argument is required in the constuctor for JopLnProp3DAcoIsoDenQ_DEO2_FDTD but not JopNlProp3DAcoIsoDenQ_DEO2_FDTD. This constuctor is shown in the last example above.
  • dtmod The sample rate for the modeled data. You can establish a lower bound for the modeling sample rate with the following expression: dt = 0.75 * 0.38 * max(dx,dz) / maximum(v). Note that the usual Courant–Friedrichs–Lewy condition (CFL condition) for stability in finite difference modeling is modified heuristically to include the impact of the visco-acoustic implementation of this operator, requiring a 25% smaller smaple rate.
  • dtrec The number of time samples in the recorded data. Note requirement that dtrec be an even multiple of dtmod.
  • ntrec The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined from ntrec, dtrec, and dtmod.

Optional Parameters

Defaults for arguments are shown inside square brackets.

  • b [ones(Float32,0,0,0)] the buoyancy (reciprocal density) array.
  • srcfieldfile [joinpath(tempdir(), "field-977465bb-cd8e-46ef-91e2-49b567797c0b.bin")] the full path to a scratch file used for the serializationof the compressed nonlinear source wavefield.
  • comptype [nothing] the type of compression to use for the serialization of the nonlinear source wavefield. The type of compression must be one of:
    • nothing - no compression.
    • Float32 - if wavefield is Float64, then do a simple conversion to Float32.
    • UInt32 - compression using CvxCompress (windowing + 2D wavelet transform + thresholding + quantization + run-length-encoding).
  • compscale [1e-2] determines the thresholding for the compression of the nonlinear source wavefield prior to serialization. Larger values mean more aggressive compression. You can likely increase compscale to 1.0 before you start to notice differences in output from Jacobian operations.
  • nz_subcube, ny_subcube, nx_subcube [32] The Z, Y, and X sizes of windows used for compression of the nonlinear source wavefield with CvxCompress. Note the requirement [8 <= n*_subcube <=256].
  • isinterior [false] boolean flag that indicates how the nonlinear source wavefield is saved. For large models, operation will be faster with isinterior = true, but the linearization correctness test may fail.
    • true the entire model including absorbing boundaries is serialized and deserialized
    • false the interior part of the model excluding absorbing boundaries is serialized and deserialized
  • sz [0.0] Array of source Z coordinate. Note that if multiple sources are provided, they will will be injected simultaneously during finite difference evolution.
  • sy [0.0] Array of source Y coordinate.
  • sx [0.0] Array of source X coordinate.
  • st [0.0] Array of source delay times.
  • interpmethod [:hicks] Type of physical interpolation for sources and receivers. For locations that are not on the physical grid coordinates, interpolation is used either to inject or extract information. interpmethod must be one of:
    • :hicks Hicks 3D sinc interpolation (up to 8x8x8 nonzero points per location)
    • :linear bilinear interpolation (up to 2x2x2 nonzero points per location)
  • rz [[0.0]] 2D array of receiver Z coordinates
  • ry [[0.0]] 2D array of receiver Y coordinates
  • rx [[0.0]] 2D array of receiver Z coordinates
  • z0 [0.0] Origin of physical coordinates for the Z dimension.
  • y0 [0.0] Origin of physical coordinates for the Y dimension.
  • x0 [0.0] Origin of physical coordinates for the X dimension.
  • dz [10.0] Spacing of physical coordinates in the Z dimension.
  • dy [10.0] Spacing of physical coordinates in the Y dimension.
  • dx [10.0] Spacing of physical coordinates in the X dimension.
  • freqQ [5.0] The center frequency for the Maxwell body approximation to dissipation only attenuation. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qMin [0.1] The minimum value for Qp at the boundary of the model used in our Maxwell body approximation to dissipation only attenuation. This is not a physically meaningful value for Qp, as we use the attenuation to implement absorbing boundary conditions and eliminate outgoing waves on the boundaries of the computational domain. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qInterior [100.0] the value for Qp in the interior of the model used in our Maxwell body approximation to dissipation only attenuation. This is the value for Qp away from the absorbing boundaries and is a physically meaningful value. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • padz, pady, padx [0.0], [0.0], [0.0] - apply extra padding to the survey determined aperture in Ginsu. Please see Ginsu for more information.
  • nbz_cache, nby_cache, nbx_cache [512], [8], [8] The size of cache blocks in the Z, X, and Y dimensions. In general the cache block in the Z (fast) dimension should be ≥ the entire size of that dimension, and the cache block size in the slower dimensions is generally small in order to allow the entire block to fit in cache.
  • nsponge [50] The number of grid cells to use for the absorbing boundary. For high fidelity modeling this should be > 60 grid points, but can be significantly smaller for some use cases like low frequency full waveform inversion.
  • wavelet [WaveletCausalRicker(f=5.0)] The source wavelet, can be specified as either a Wavelet type or an array.
  • freesurface [false] Determines if a free surface (true) or absorbing (false) top boundary condition is applied.
  • imgcondition ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint only exists for "standard" imaging condition currently.
  • RTM_weight [0.5] determines the balance of short wavelengths and long wavelengths in the imaging condition. A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition.
  • nthreads [Sys.CPU_THREADS] The number of threads to use for OpenMP parallelization of the modeling.
  • reportinterval [500] The interval at which information about the propagtion is logged.

See also: Ginsu, WaveletSine, WaveletRicker, WaveletMinPhaseRicker, WaveletDerivRicker, WaveletCausalRicker, WaveletOrmsby, WaveletMinPhaseOrmsby

source
JetPackWaveFD.JopLnProp3DAcoTTIDenQ_DEO2_FDTD — Method
JopNlProp3DAcoTTIDenQ_DEO2_FDTD(; kwargs...)
JopLnProp3DAcoTTIDenQ_DEO2_FDTD(; m₀, kwargs...)

Create a Jets nonlinear or linearized operator for 3D pseudo- visco-acoustic, tilted transverse isotropy, variable density modeling.

Model Parameters

This propagator operates with four model parameters, as shown in the table below.

ParameterDescription
vP wave velocity
ϵModified Thomsen's weak anisotropy parameter
ηModified Alkhalifah's weak anisotropy parameter η
bbuoyancy (reciprocal density)
θsymmetry axis tilt angle from vertical (radians)
ϕsymmetry axis aziumuth angle CCW from x (radians)

Pseudo-acoustic approximation

With the pseudo-acoustic approximation we employ, the transformation to self-adjoint form requires a parameter f representing the average ratio of shear-wave to p-wave velocity, and a modfication to Alkhalifah's weak anisotropy parameter η. We show formulas for these two parameters below. The parameter f is specified as a scalar for the entire model, with default value 0.85, implying a shear velocity around 38% of the p wave velocity.

  • η = sqrt[2(ϵ-δ)/(f+2ϵ)]
  • f = 1 - Vₛ² / Vₚ²

For more information about the pseudo-acoustic approximation we employ please see https://library.seg.org/doi/10.1190/segam2016-13878451.1.

Active and passive model parameters

There are two modes of operation that define different sets of active and passive parameters, velocity only and velocity and anisotropy.

An active parameter can be inverted for using the Jacobian linearized operator machinery, and a passive parameter is constant. The two modes are:

ModeActive ParametersPassive Parameters
velocity only[v][ϵ, η, b, θ, ϕ ]
velocity and anisotropy[v, ϵ, η][b, θ, ϕ ]

To make a parameter passive, you pass the array for that parameter in the constructor for the operator, implying that it is in state and does not change with the action of the operator. Parameters that are not passed to the constructor are assumed to be active, and will be part of the model that the operator acts on, stored as a 3D array with the following indices: [Z coord][X coord][Active Parameter]

Examples

Model and acquisition geometry setup

  1. load modules Jets, WaveFD, and JetPackWaveFDFDFD
  2. set up the model discretization, coordinate size and spacing
  3. set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
  4. create constant models for v, ϵ, η, b, θ, ϕ
using Jets, WaveFD, JetPackWaveFDFDFD
nz,ny,nx = 100, 80, 60                        # spatial discretization size
dz,dy,dx = 20.0, 20.0, 20.0                   # spatial discretization sampling
ntrec = 1101                                  # number of temporal samples in recorded data
dtrec = 0.0100                                # sample rate for recorded data
dtmod = 0.0025                                # sample rate for modeled data
wavelet = WaveletCausalRicker(f=5.0)          # type of wavelet to use (source signature) 
sz = dz                                       # source Z location 
sy = dy*(ny/2)                                # source Y location 
sx = dx*(nx/2)                                # source X location 
rz = [dz for iy = 1:ny, ix = 1:nx][:];        # Array of receiver Y locations 
ry = [(iy-1)*dy for iy = 1:ny, ix = 1:nx][:]; # Array of receiver Y locations
rx = [(ix-1)*dx for iy = 1:ny, ix=1:nx][:];   # Array of receiver X locations
b = ones(Float32, nz, ny, nx);                # constant buoyancy
ϵ = 0.1 * ones(Float32, nz, ny, nx);          # constant epsilon
η = 0.2 * ones(Float32, nz, ny, nx);          # constant eta
θ = (π/8) * ones(Float32, nz, ny, nx);        # constant tilt angle
ϕ = (π/3) * ones(Float32, nz, ny, nx);        # constant azimuth angle

Construct and apply the nonlinear operator (v is active; ϵ, η, b, θ, ϕ are passive)

  1. create the nonlinear operator F
  2. create the constant velocity model m₀
  3. perform nonlinear forward modeling with constant velocity model m₀ and return the resulting modeled data in d
F = JopNlProp3DAcoTTIDenQ_DEO2_FDTD(; b = b, ϵ = ϵ, η = η, θ = θ, ϕ = ϕ, isinterior=true, nsponge = 10, 
    ntrec = ntrec, dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, 
    wavelet = wavelet, sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,:,1] .= 1500;
d = F*m₀;              # forward nonlinear op

Construct and apply the nonlinear operator (v, ϵ, η are active; b, θ, ϕ are passive)

  1. create the nonlinear operator F
  2. create the model vector and set parameters v, ϵ, η
  3. perform nonlinear forward modeling with model m₀ and return the resulting modeled data in d
F = JopNlProp3DAcoTTIDenQ_DEO2_FDTD(; b = b, θ = θ, ϕ = ϕ, isinterior=true, nsponge = 10, 
    ntrec = ntrec, dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, 
    wavelet = wavelet, sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,:,1] .= 1500;
m₀[:,:,:,2] .= 0.1;
m₀[:,:,:,3] .= 0.2;
d = F*m₀;              # forward nonlinear op

Construct and apply the linearized Jacobian operator (method 1)

For this example we assume in the model that v, ϵ, η are active parameters, and b, θ, ϕ are passive.

  1. create the nonlinear operator F directly by construction of a jet.
  2. create the model vector and set parameters v, ϵ, η
  3. create the Jacobian operator J by linearization of F at point m₀
  4. create a random model perturbation vector δm.
  5. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  6. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm.

Note that the Jacobian operators J and J' require the serialized nonlinear forward wavefield, and this is generated automatically whenever required. If you watch the logging to standard out from this example, you will first see the finite difference evolution for the nonlinear forward, followed by the linearized forward, and finally the linearized adjoint.

F = JopNlProp3DAcoTTIDenQ_DEO2_FDTD(; b = b, θ = θ, ϕ = ϕ, isinterior=true, nsponge = 10, 
    ntrec = ntrec, dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,:,1] .= 1500;
m₀[:,:,:,2] .= 0.1;
m₀[:,:,:,3] .= 0.2;
J = jacobian(F, m₀)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Construct and apply the linearized Jacobian operator (method 2)

For this example we assume in the model that v, ϵ, η are active parameters, and b, θ, ϕ are passive.

  1. create the constant velocity model m₀
  2. create the Jacobian operator J at point m₀ directly by construction of a jet.
  3. create a random model perturbation vector δm.
  4. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  5. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm. 
m₀ = zeros(Float32, nz, ny, nx, 3);
m₀[:,:,:,1] .= 1500;
m₀[:,:,:,2] .= 0.1;
m₀[:,:,:,3] .= 0.2;
J = JopLnProp3DAcoTTIDenQ_DEO2_FDTD(; m₀ = m₀, b = b, θ = θ, ϕ = ϕ, isinterior=true, nsponge = 10, 
    ntrec = ntrec, dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, 
    wavelet = wavelet, sz = sz, sy = sy, sx = sx, ry = ry, rz = rz, rx = rx)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Required Parameters

  • m₀ the point at which the jet is linearized. Note this argument is required in the constuctor for JopLnProp2DAcoTTIDenQ_DEO2_FDTD but not JopNlProp2DAcoTTIDenQ_DEO2_FDTD. This constuctor is shown in the last example above. Please note that you must consider which parameters are active and passive, per the discussion on model parameters above and examples.
  • dtmod The sample rate for the modeled data. You can establish a lower bound for the modeling sample rate with the following expression: dt = 0.75 * 0.38 * max(dx,dz) / maximum(v). Note that the usual Courant–Friedrichs–Lewy condition (CFL condition) for stability in finite difference modeling is modified heuristically to include the impact of the visco-acoustic implementation of this operator, requiring a 25% smaller smaple rate.
  • dtrec The number of time samples in the recorded data. Note requirement that dtrec be an even multiple of dtmod.
  • ntrec The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined from ntrec, dtrec, and dtmod.

Optional Parameters

Defaults for arguments are shown inside square brackets.

  • b [ones(Float32,0,0,0)] the buoyancy (reciprocal density) array.
  • ϕ [zeros(Float32,0,0,0)] the symmetry axies azimuth angle counter-clockwise from the x axis, in radians. Assumed to be zero if not specified.
  • θ [zeros(Float32,0,0,0)] the symmetry axies tilt angle in radians. Assumed to be zero if not specified. θ must be an array the same size as buoyancy b.
  • f [0.85] See the discussion above regarding the Pseudo-acoustic approximation used here
  • srcfieldfile ["field-4fee5977-1631-499f-9f65-93b02219bb06.bin"] the full path to a scratch file used for the serializationof the compressed nonlinear source wavefield.
  • comptype [nothing] the type of compression to use for the serialization of the nonlinear source wavefield. The type of compression must be one of:
    • nothing - no compression.
    • Float32 - if wavefield is Float64, then do a simple conversion to Float32.
    • UInt32 - compression using CvxCompress (windowing + 2D wavelet transform + thresholding + quantization + run-length-encoding).
  • compscale [1e-2] determines the thresholding for the compression of the nonlinear source wavefield prior to serialization. Larger values mean more aggressive compression. You can likely increase compscale to 1.0 before you start to notice differences in output from Jacobian operations.
  • nz_subcube, ny_subcube, nx_subcube [32] The Z, Y, and X sizes of windows used for compression of the nonlinear source wavefield with CvxCompress. Note the requirement [8 <= n*_subcube <=256].
  • isinterior [false] boolean flag that indicates how the nonlinear source wavefield is saved. For large models, operation will be faster with isinterior = true, but the linearization correctness test may fail.
    • true the entire model including absorbing boundaries is serialized and deserialized
    • false the interior part of the model excluding absorbing boundaries is serialized and deserialized
  • sz [0.0] Array of source Z coordinate. Note that if multiple sources are provided, they will will be injected simultaneously during finite difference evolution.
  • sy [0.0] Array of source Y coordinate.
  • sx [0.0] Array of source X coordinate.
  • st [0.0] Array of source delay times.
  • interpmethod [:hicks] Type of physical interpolation for sources and receivers. For locations that are not on the physical grid coordinates, interpolation is used either to inject or extract information. interpmethod must be one of:
    • :hicks Hicks 3D sinc interpolation (up to 8x8x8 nonzero points per location)
    • :linear bilinear interpolation (up to 2x2x2 nonzero points per location)
  • rz [[0.0]] 2D array of receiver Z coordinates
  • ry [[0.0]] 2D array of receiver Y coordinates
  • rx [[0.0]] 2D array of receiver Z coordinates
  • z0 [0.0] Origin of physical coordinates for the Z dimension.
  • y0 [0.0] Origin of physical coordinates for the Y dimension.
  • x0 [0.0] Origin of physical coordinates for the X dimension.
  • dz [10.0] Spacing of physical coordinates in the Z dimension.
  • dy [10.0] Spacing of physical coordinates in the Y dimension.
  • dx [10.0] Spacing of physical coordinates in the X dimension.
  • freqQ [5.0] The center frequency for the Maxwell body approximation to dissipation only attenuation. Please see JetPackWaveFDFDFD package documentation for more information concerning the attenuation model.
  • qMin [0.1] The minimum value for Qp at the boundary of the model used in our Maxwell body approximation to dissipation only attenuation. This is not a physically meaningful value for Qp, as we use the attenuation to implement absorbing boundary conditions and eliminate outgoing waves on the boundaries of the computational domain. Please see JetPackWaveFDFDFD package documentation for more information concerning the attenuation model.
  • qInterior [100.0] the value for Qp in the interior of the model used in our Maxwell body approximation to dissipation only attenuation. This is the value for Qp away from the absorbing boundaries and is a physically meaningful value. Please see JetPackWaveFDFDFD package documentation for more information concerning the attenuation model.
  • padz, pady, padx [0.0], [0.0], [0.0] - apply extra padding to the survey determined aperture in Ginsu. Please see Ginsu for more information.
  • nbz_cache, nby_cache, nbx_cache [512], [8], [8] The size of cache blocks in the Z, X, and Y dimensions. In general the cache block in the Z (fast) dimension should be ≥ the entire size of that dimension, and the cache block size in the slower dimensions is generally small in order to allow the entire block to fit in cache.
  • nsponge [50] The number of grid cells to use for the absorbing boundary. For high fidelity modeling this should be > 60 grid points, but can be significantly smaller for some use cases like low frequency full waveform inversion.
  • wavelet [WaveletCausalRicker(f=5.0)] The source wavelet, can be specified as either a Wavelet type or an array.
  • freesurface [false] Determines if a free surface (true) or absorbing (false) top boundary condition is applied.
  • imgcondition ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint only exists for "standard" imaging condition currently.
  • RTM_weight [0.5] determines the balance of short wavelengths and long wavelengths in the imaging condition. A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition.
  • nthreads [Sys.CPU_THREADS] The number of threads to use for OpenMP parallelization of the modeling.
  • reportinterval [500] The interval at which information about the propagtion is logged.

See also: Ginsu, WaveletSine, WaveletRicker, WaveletMinPhaseRicker, WaveletDerivRicker, WaveletCausalRicker, WaveletOrmsby, WaveletMinPhaseOrmsby

source
JetPackWaveFD.JopLnProp3DAcoVTIDenQ_DEO2_FDTD — Method
JopNlProp3DAcoVTIDenQ_DEO2_FDTD(; kwargs...)
JopLnProp3DAcoVTIDenQ_DEO2_FDTD(; m₀, kwargs...)

Create a Jets nonlinear or linearized operator for 3D pseudo- visco-acoustic, vertical transverse isotropy, variable density modeling.

Model Parameters

This propagator operates with four model parameters, as shown in the table below.

ParameterDescription
vP wave velocity
ϵModified Thomsen's weak anisotropy parameter
ηModified Alkhalifah's weak anisotropy parameter η
bbuoyancy (reciprocal density)

Pseudo-acoustic approximation

With the pseudo-acoustic approximation we employ, the transformation to self-adjoint form requires a parameter f representing the average ratio of shear-wave to p-wave velocity, and a modfication to Alkhalifah's weak anisotropy parameter η. We show formulas for these two parameters below. The parameter f is specified as a scalar for the entire model, with default value 0.85, implying a shear velocity around 38% of the p wave velocity.

  • η = sqrt[2(ϵ-δ)/(f+2ϵ)]
  • f = 1 - Vₛ² / Vₚ²

For more information about the pseudo-acoustic approximation we employ please see https://library.seg.org/doi/10.1190/segam2016-13878451.1.

Active and passive model parameters

There are two modes of operation that define different sets of active and passive parameters, velocity only and velocity and anisotropy.

An active parameter can be inverted for using the Jacobian linearized operator machinery, and a passive parameter is constant. The two modes are:

ModeActive ParametersPassive Parameters
velocity only[v][ϵ, η, b]
velocity and anisotropy[v, ϵ, η][b]

To make a parameter passive, you pass the array for that parameter in the constructor for the operator, implying that it is in state and does not change with the action of the operator. Parameters that are not passed to the constructor are assumed to be active, and will be part of the model that the operator acts on, stored as a 3D array with the following indices: [Z coord][X coord][Active Parameter]

Examples

Model and acquisition geometry setup

  1. load modules Jets, WaveFD, and JetPackWaveFD
  2. set up the model discretization, coordinate size and spacing
  3. set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
  4. create constant models for v, ϵ, η, b
using Jets, WaveFD, JetPackWaveFD
nz,ny,nx = 100, 80, 60                        # spatial discretization size
dz,dy,dx = 20.0, 20.0, 20.0                   # spatial discretization sampling
ntrec = 1101                                  # number of temporal samples in recorded data
dtrec = 0.0100                                # sample rate for recorded data
dtmod = 0.0025                                # sample rate for modeled data
wavelet = WaveletCausalRicker(f=5.0)          # type of wavelet to use (source signature) 
sz = dz                                       # source Z location 
sy = dy*(ny/2)                                # source Y location 
sx = dx*(nx/2)                                # source X location 
rz = [dz for iy = 1:ny, ix = 1:nx][:];        # Array of receiver Y locations 
ry = [(iy-1)*dy for iy = 1:ny, ix = 1:nx][:]; # Array of receiver Y locations
rx = [(ix-1)*dx for iy = 1:ny, ix=1:nx][:];   # Array of receiver X locations
b = ones(Float32, nz, ny, nx);                # constant buoyancy
ϵ = 0.1 * ones(Float32, nz, ny, nx);          # constant epsilon
η = 0.2 * ones(Float32, nz, ny, nx);          # constant eta

Construct and apply the nonlinear operator (v is active; ϵ, η, b are passive)

  1. create the nonlinear operator F
  2. create the constant velocity model m₀
  3. perform nonlinear forward modeling with constant velocity model m₀ and return the resulting modeled data in d
F = JopNlProp3DAcoVTIDenQ_DEO2_FDTD(; b = b, ϵ = ϵ, η = η, isinterior=true, nsponge = 10, 
    ntrec = ntrec, dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, 
    wavelet = wavelet, sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,:,1] .= 1500;
d = F*m₀;              # forward nonlinear op

Construct and apply the nonlinear operator (v, ϵ, η are active; b is passive)

  1. create the nonlinear operator F
  2. create the model vector and set parameters v, ϵ, η
  3. perform nonlinear forward modeling with model m₀ and return the resulting modeled data in d
F = JopNlProp3DAcoVTIDenQ_DEO2_FDTD(; b = b, isinterior=true, nsponge = 10, 
    ntrec = ntrec, dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, 
    wavelet = wavelet, sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,:,1] .= 1500;
m₀[:,:,:,2] .= 0.1;
m₀[:,:,:,3] .= 0.2;
d = F*m₀;              # forward nonlinear op

Construct and apply the linearized Jacobian operator (method 1)

For this example we assume in the model that v, ϵ, η are active parameters, and b is passive.

  1. create the nonlinear operator F directly by construction of a jet.
  2. create the model vector and set parameters v, ϵ, η
  3. create the Jacobian operator J by linearization of F at point m₀
  4. create a random model perturbation vector δm.
  5. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  6. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm.

Note that the Jacobian operators J and J' require the serialized nonlinear forward wavefield, and this is generated automatically whenever required. If you watch the logging to standard out from this example, you will first see the finite difference evolution for the nonlinear forward, followed by the linearized forward, and finally the linearized adjoint.

F = JopNlProp3DAcoVTIDenQ_DEO2_FDTD(; b = b, isinterior=true, nsponge = 10, ntrec = ntrec, 
    dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,:,1] .= 1500;
m₀[:,:,:,2] .= 0.1;
m₀[:,:,:,3] .= 0.2;
J = jacobian(F, m₀)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Construct and apply the linearized Jacobian operator (method 2)

For this example we assume in the model that v, ϵ, η are active parameters, and b is passive.

  1. create the constant velocity model m₀
  2. create the Jacobian operator J at point m₀ directly by construction of a jet.
  3. create a random model perturbation vector δm.
  4. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  5. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm. 
m₀ = zeros(Float32, nz, ny, nx, 3);
m₀[:,:,:,1] .= 1500;
m₀[:,:,:,2] .= 0.1;
m₀[:,:,:,3] .= 0.2;
J = JopLnProp3DAcoVTIDenQ_DEO2_FDTD(; m₀ = m₀, b = b, isinterior=true, nsponge = 10, 
    ntrec = ntrec, dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, 
    wavelet = wavelet, sz = sz, sy = sy, sx = sx, ry = ry, rz = rz, rx = rx)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Required Parameters

  • m₀ the point at which the jet is linearized. Note this argument is required in the constuctor for JopLnProp2DAcoVTIDenQ_DEO2_FDTD but not JopNlProp2DAcoVTIDenQ_DEO2_FDTD. This constuctor is shown in the last example above. Please note that you must consider which parameters are active and passive, per the discussion on model parameters above and examples.
  • dtmod The sample rate for the modeled data. You can establish a lower bound for the modeling sample rate with the following expression: dt = 0.75 * 0.38 * max(dx,dz) / maximum(v). Note that the usual Courant–Friedrichs–Lewy condition (CFL condition) for stability in finite difference modeling is modified heuristically to include the impact of the visco-acoustic implementation of this operator, requiring a 25% smaller smaple rate.
  • dtrec The number of time samples in the recorded data. Note requirement that dtrec be an even multiple of dtmod.
  • ntrec The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined from ntrec, dtrec, and dtmod.

Optional Parameters

Defaults for arguments are shown inside square brackets.

  • b [ones(Float32,0,0,0)] the buoyancy (reciprocal density) array.
  • f [0.85] See the discussion above regarding the Pseudo-acoustic approximation used here
  • srcfieldfile [joinpath(tempdir(), "field-e9ebc6be-af0c-4c5a-b7dc-5b5ccb6774a4.bin")] the full path to a scratch file used for the serializationof the compressed nonlinear source wavefield.
  • comptype [nothing] the type of compression to use for the serialization of the nonlinear source wavefield. The type of compression must be one of:
    • nothing - no compression.
    • Float32 - if wavefield is Float64, then do a simple conversion to Float32.
    • UInt32 - compression using CvxCompress (windowing + 2D wavelet transform + thresholding + quantization + run-length-encoding).
  • compscale [1e-2] determines the thresholding for the compression of the nonlinear source wavefield prior to serialization. Larger values mean more aggressive compression. You can likely increase compscale to 1.0 before you start to notice differences in output from Jacobian operations.
  • nz_subcube, ny_subcube, nx_subcube [32] The Z, Y, and X sizes of windows used for compression of the nonlinear source wavefield with CvxCompress. Note the requirement [8 <= n*_subcube <=256].
  • isinterior [false] boolean flag that indicates how the nonlinear source wavefield is saved. For large models, operation will be faster with isinterior = true, but the linearization correctness test may fail.
    • true the entire model including absorbing boundaries is serialized and deserialized
    • false the interior part of the model excluding absorbing boundaries is serialized and deserialized
  • sz [0.0] Array of source Z coordinate. Note that if multiple sources are provided, they will will be injected simultaneously during finite difference evolution.
  • sy [0.0] Array of source Y coordinate.
  • sx [0.0] Array of source X coordinate.
  • st [0.0] Array of source delay times.
  • interpmethod [:hicks] Type of physical interpolation for sources and receivers. For locations that are not on the physical grid coordinates, interpolation is used either to inject or extract information. interpmethod must be one of:
    • :hicks Hicks 3D sinc interpolation (up to 8x8x8 nonzero points per location)
    • :linear bilinear interpolation (up to 2x2x2 nonzero points per location)
  • rz [[0.0]] 2D array of receiver Z coordinates
  • ry [[0.0]] 2D array of receiver Y coordinates
  • rx [[0.0]] 2D array of receiver Z coordinates
  • z0 [0.0] Origin of physical coordinates for the Z dimension.
  • y0 [0.0] Origin of physical coordinates for the Y dimension.
  • x0 [0.0] Origin of physical coordinates for the X dimension.
  • dz [10.0] Spacing of physical coordinates in the Z dimension.
  • dy [10.0] Spacing of physical coordinates in the Y dimension.
  • dx [10.0] Spacing of physical coordinates in the X dimension.
  • freqQ [5.0] The center frequency for the Maxwell body approximation to dissipation only attenuation. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qMin [0.1] The minimum value for Qp at the boundary of the model used in our Maxwell body approximation to dissipation only attenuation. This is not a physically meaningful value for Qp, as we use the attenuation to implement absorbing boundary conditions and eliminate outgoing waves on the boundaries of the computational domain. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qInterior [100.0] the value for Qp in the interior of the model used in our Maxwell body approximation to dissipation only attenuation. This is the value for Qp away from the absorbing boundaries and is a physically meaningful value. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • padz, pady, padx [0.0], [0.0], [0.0] - apply extra padding to the survey determined aperture in Ginsu. Please see Ginsu for more information.
  • nbz_cache, nby_cache, nbx_cache [512], [8], [8] The size of cache blocks in the Z, X, and Y dimensions. In general the cache block in the Z (fast) dimension should be ≥ the entire size of that dimension, and the cache block size in the slower dimensions is generally small in order to allow the entire block to fit in cache.
  • nsponge [50] The number of grid cells to use for the absorbing boundary. For high fidelity modeling this should be > 60 grid points, but can be significantly smaller for some use cases like low frequency full waveform inversion.
  • wavelet [WaveletCausalRicker(f=5.0)] The source wavelet, can be specified as either a Wavelet type or an array.
  • freesurface [false] Determines if a free surface (true) or absorbing (false) top boundary condition is applied.
  • imgcondition ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint only exists for "standard" imaging condition currently.
  • RTM_weight [0.5] determines the balance of short wavelengths and long wavelengths in the imaging condition. A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition.
  • nthreads [Sys.CPU_THREADS] The number of threads to use for OpenMP parallelization of the modeling.
  • reportinterval [500] The interval at which information about the propagtion is logged.

See also: Ginsu, WaveletSine, WaveletRicker, WaveletMinPhaseRicker, WaveletDerivRicker, WaveletCausalRicker, WaveletOrmsby, WaveletMinPhaseOrmsby

source
JetPackWaveFD.JopNlProp2DAcoIsoDenQ_DEO2_FDTD — Method
JopNlProp2DAcoIsoDenQ_DEO2_FDTD(; kwargs...)
JopLnProp2DAcoIsoDenQ_DEO2_FDTD(; v, kwargs...)

Create a Jets nonlinear or linearized operator for 2D visco-acoustic, isotropic, variable density modeling.

Model Parameters

This propagator operates with two model parameters, as shown in the table below.

ParameterDescription
vP wave velocity
bbuoyancy (reciprocal density)

Velocity v and buoyancy b can be active parameters and inverted for using the Jacobian machinery, or passive parameters that are constant. If a parameter is active there will be a component of model assigned to it.

The model is 3D with the slowest dimension for model parameter. The earth model property maps to the slow- dimension array index via the state(F).active_modelset dictionary.

Examples

Model and acquisition geometry setup

  1. load modules Jets, WaveFD, and JetPackWaveFD
  2. set up the model discretization, coordinate size and spacing
  3. set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
  4. create constant buoyancy array
using Jets, WaveFD, JetPackWaveFD
nz,nx = 200, 100                      # spatial discretization size
dz,dx = 20.0, 20.0                    # spatial discretization sampling
ntrec = 1101                          # number of temporal samples in recorded data
dtrec = 0.0100                        # sample rate for recorded data
dtmod = 0.0025                        # sample rate for modeled data
wavelet = WaveletCausalRicker(f=5.0)  # type of wavelet to use (source signature) 
sz = dz                               # source Z location 
sx = dx*(nx/2)                        # source X location 
rx = dx*[0:0.5:nx-1;];                # array of receiver X locations
rz = 2*dz*ones(length(0:0.5:nx-1));   # array of receiver Z locations
b = ones(Float32, nz, nx);            # buoyancy model (reciprocal density)

Construct and apply the nonlinear operator

  1. create the nonlinear operator F
  2. create the constant velocity model m₀
  3. perform nonlinear forward modeling with constant velocity model m₀ and return the resulting modeled data in d
F = JopNlProp2DAcoIsoDenQ_DEO2_FDTD(; b = b, ntrec = ntrec, 
    dz = dz, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = 1500 .* ones(domain(F));
d = F*m₀;              # forward nonlinear op

Construct and apply the linearized Jacobian operator (method 1)

  1. create the nonlinear operator F directly by construction of a jet.
  2. create the constant velocity model m₀
  3. create the Jacobian operator J by linearization of F at point m₀
  4. create a random model perturbation vector δm.
  5. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  6. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm.

Note that the Jacobian operators J and J' require the serialized nonlinear forward wavefield, and this is generated automatically whenever required. If you watch the logging to standard out from this example, you will first see the finite difference evolution for the nonlinear forward, followed by the linearized forward, and finally the linearized adjoint.

F = JopNlProp2DAcoIsoDenQ_DEO2_FDTD(; b = b, ntrec = ntrec, 
    dz = dz, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = 1500 .* ones(domain(F));
J = jacobian(F, m₀)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Construct and apply the linearized Jacobian operator (method 2)

  1. create the constant velocity model m₀
  2. create the Jacobian operator J at point m₀ directly by construction of a jet.
  3. create a random model perturbation vector δm.
  4. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  5. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm. 
m₀ = 1500 * ones(Float32, nz, nx);
J = JopLnProp2DAcoIsoDenQ_DEO2_FDTD(; v = m₀, b = b, ntrec = ntrec, 
    dz = dz, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sx = sx, rz = rz, rx = rx)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Required Parameters

  • v the point at which the jet is linearized. Note this argument is required in the constuctor for JopLnProp2DAcoIsoDenQ_DEO2_FDTD but not JopNlProp2DAcoIsoDenQ_DEO2_FDTD. This constuctor is shown in the last example above.
  • dtmod The sample rate for the modeled data. You can establish a lower bound for the modeling sample rate with the following expression: dt = 0.75 * 0.38 * max(dx,dz) / maximum(v). Note that the usual Courant–Friedrichs–Lewy condition (CFL condition) for stability in finite difference modeling is modified heuristically to include the impact of the visco-acoustic implementation of this operator, requiring a 25% smaller smaple rate.
  • dtrec The number of time samples in the recorded data. Note requirement that dtrec be an even multiple of dtmod.
  • ntrec The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined from ntrec, dtrec, and dtmod.

Optional Parameters

Defaults for arguments are shown inside square brackets.

  • b [ones(Float32,0,0)] the buoyancy (reciprocal density) array.
  • srcfieldfile [joinpath(tempdir(), "field-c3628917-2e59-43dc-9fb7-c58a2f8e1d4d.bin")] the full path to a scratch file used for the serializationof the compressed nonlinear source wavefield.
  • comptype [nothing] the type of compression to use for the serialization of the nonlinear source wavefield. The type of compression must be one of:
    • nothing - no compression.
    • Float32 - if wavefield is Float64, then do a simple conversion to Float32.
    • UInt32 - compression using CvxCompress (windowing + 2D wavelet transform + thresholding + quantization + run-length-encoding).
  • compscale [1e-2] determines the thresholding for the compression of the nonlinear source wavefield prior to serialization. Larger values mean more aggressive compression. You can likely increase compscale to 1.0 before you start to notice differences in output from Jacobian operations.
  • nz_subcube, nx_subcube [32] The Z and X sizes of windows used for compression of the nonlinear source wavefield with CvxCompress. Note the requirement [8 <= n*_subcube <=256].
  • isinterior [false] boolean flag that indicates how the nonlinear source wavefield is saved. For large models, operation will be faster with isinterior = true, but the linearization correctness test may fail.
    • true the entire model including absorbing boundaries is serialized and deserialized
    • false the interior part of the model excluding absorbing boundaries is serialized and deserialized
  • sz [0.0] Array of source Z coordinate. Note that if multiple sources are provided, they will will be injected simultaneously during finite difference evolution.
  • sx [0.0] Array of source X coordinate.
  • st [0.0] Array of source delay times.
  • interpmethod [:hicks] Type of physical interpolation for sources and receivers. For locations that are not on the physical grid coordinates, interpolation is used either to inject or extract information. interpmethod must be one of:
    • :hicks Hicks 2D sinc interpolation (up to 8x8 nonzero points per location)
    • :linear bilinear interpolation (up to 2x2 nonzero points per location)
  • rz [[0.0]] 2D array of receiver Z coordinates
  • rx [[0.0]] 2D array of receiver Z coordinates
  • z0 [0.0] Origin of physical coordinates for the Z dimension.
  • x0 [0.0] Origin of physical coordinates for the X dimension.
  • dz [10.0] Spacing of physical coordinates in the Z dimension.
  • dx [10.0] Spacing of physical coordinates in the X dimension.
  • freqQ [5.0] The center frequency for the Maxwell body approximation to dissipation only attenuation. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qMin [0.1] The minimum value for Qp at the boundary of the model used in our Maxwell body approximation to dissipation only attenuation. This is not a physically meaningful value for Qp, as we use the attenuation to implement absorbing boundary conditions and eliminate outgoing waves on the boundaries of the computational domain. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qInterior [100.0] the value for Qp in the interior of the model used in our Maxwell body approximation to dissipation only attenuation. This is the value for Qp away from the absorbing boundaries and is a physically meaningful value. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • padz, padx [0.0], [0.0] - apply extra padding to the survey determined aperture in Ginsu. Please see Ginsu for more information.
  • nbz_cache, nbx_cache [512], [8] The size of cache blocks in the Z and X dimensions. In general the cache block in the Z (fast) dimension should be ≥ the entire size of that dimension, and the cache block size in the slower dimension is generally small in order to allow the entire block to fit in cache.
  • nsponge [50] The number of grid cells to use for the absorbing boundary. For high fidelity modeling this should be > 60 grid points, but can be significantly smaller for some use cases like low frequency full waveform inversion.
  • wavelet [WaveletCausalRicker(f=5.0)] The source wavelet, can be specified as either a Wavelet type or an array.
  • freesurface [false] Determines if a free surface (true) or absorbing (false) top boundary condition is applied.
  • imgcondition ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint only exists for "standard" imaging condition currently.
  • RTM_weight [0.5] determines the balance of short wavelengths and long wavelengths in the imaging condition. A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition.
  • nthreads [Sys.CPU_THREADS] The number of threads to use for OpenMP parallelization of the modeling.
  • reportinterval [500] The interval at which information about the propagtion is logged.

See also: Ginsu, WaveletSine, WaveletRicker, WaveletMinPhaseRicker, WaveletDerivRicker, WaveletCausalRicker, WaveletOrmsby, WaveletMinPhaseOrmsby

source
JetPackWaveFD.JopNlProp2DAcoTTIDenQ_DEO2_FDTD — Method
JopNlProp2DAcoTTIDenQ_DEO2_FDTD(; kwargs...)
JopLnProp2DAcoTTIDenQ_DEO2_FDTD(; m₀, kwargs...)

Create a Jets nonlinear or linearized operator for 2D pseudo- visco-acoustic, tilted transverse isotropy, variable density modeling.

Model Parameters

This propagator operates with five model parameters, as shown in the table below.

ParameterDescription
vP wave velocity
ϵModified Thomsen's weak anisotropy parameter
ηModified Alkhalifah's weak anisotropy parameter η
bbuoyancy (reciprocal density)
θsymmetry axis tilt angle from vertical (radians)

Pseudo-acoustic approximation

With the pseudo-acoustic approximation we employ, the transformation to self-adjoint form requires a parameter f representing the average ratio of shear-wave to p-wave velocity, and a modfication to Alkhalifah's weak anisotropy parameter η. We show formulas for these two parameters below. The parameter f is specified as a scalar for the entire model, with default value 0.85, implying a shear velocity around 38% of the p wave velocity.

  • η = sqrt[2(ϵ-δ)/(f+2ϵ)]
  • f = 1 - Vₛ² / Vₚ²

For more information about the pseudo-acoustic approximation we employ please see https://library.seg.org/doi/10.1190/segam2016-13878451.1.

Active and passive model parameters

There are two modes of operation that define different sets of active and passive parameters, velocity only and velocity and anisotropy.

An active parameter can be inverted for using the Jacobian linearized operator machinery, and a passive parameter is constant. The two modes are:

ModeActive ParametersPassive Parameters
velocity only[v][ϵ, η, b, θ ]
velocity and anisotropy[v, ϵ, η][b, θ ]

To make a parameter passive, you pass the array for that parameter in the constructor for the operator, implying that it is in state and does not change with the action of the operator. Parameters that are not passed to the constructor are assumed to be active, and will be part of the model that the operator acts on, stored as a 3D array with the following indices: [Z coord][X coord][Active Parameter]

Examples

Model and acquisition geometry setup

  1. load modules Jets, WaveFD, and JetPackWaveFD
  2. set up the model discretization, coordinate size and spacing
  3. set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
  4. create constant models for v, ϵ, η, b, θ
using Jets, WaveFD, JetPackWaveFD
nz,nx = 200, 100                      # spatial discretization size
dz,dx = 20.0, 20.0                    # spatial discretization sampling
ntrec = 1101                          # number of temporal samples in recorded data
dtrec = 0.0100                        # sample rate for recorded data
dtmod = 0.0025                        # sample rate for modeled data
wavelet = WaveletCausalRicker(f=5.0)  # type of wavelet to use (source signature) 
sz = dz                               # source Z location 
sx = dx*(nx/2)                        # source X location 
rx = dx*[0:0.5:nx-1;];                # array of receiver X locations
rz = 2*dz*ones(length(0:0.5:nx-1));   # array of receiver Z locations
b = ones(Float32, nz, nx);            # constant buoyancy
ϵ = 0.1 * ones(Float32, nz, nx);      # constant epsilon
η = 0.2 * ones(Float32, nz, nx);      # constant eta
θ = (π/8) * ones(Float32, nz, nx);    # constant tilt angle

Construct and apply the nonlinear operator (v is active; ϵ, η, b, θ are passive)

  1. create the nonlinear operator F
  2. create the model vector and set parameter v
  3. perform nonlinear forward modeling with model m₀ and return the resulting modeled data in d
F = JopNlProp2DAcoTTIDenQ_DEO2_FDTD(; b = b, ϵ = ϵ, η = η, θ = θ, ntrec = ntrec, dtrec = dtrec, 
    dtmod = dtmod, dz = dz, dx = dx, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,1] .= 1500;
d = F*m₀;              # forward nonlinear op

Construct and apply the nonlinear operator (v, ϵ, η are active; b, θ are passive)

  1. create the nonlinear operator F
  2. create the model vector and set parameters v, ϵ, η
  3. perform nonlinear forward modeling with model m₀ and return the resulting modeled data in d
F = JopNlProp2DAcoTTIDenQ_DEO2_FDTD(; b = b, θ = θ, ntrec = ntrec, dtrec = dtrec, dtmod = dtmod, 
    dz = dz, dx = dx, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,1] .= 1500;
m₀[:,:,2] .= 0.1;
m₀[:,:,3] .= 0.2;
d = F*m₀;              # forward nonlinear op

Construct and apply the linearized Jacobian operator (method 1)

For this example we assume in the model that v, ϵ, η are active parameters, and b, θ are passive.

  1. create the nonlinear operator F directly by construction of a jet.
  2. create the model vector and set parameters v, ϵ, η
  3. create the Jacobian operator J by linearization of F at point m₀
  4. create a random model perturbation vector δm.
  5. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  6. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm.

Note that the Jacobian operators J and J' require the serialized nonlinear forward wavefield, and this is generated automatically whenever required. If you watch the logging to standard out from this example, you will first see the finite difference evolution for the nonlinear forward, followed by the linearized forward, and finally the linearized adjoint.

F = JopNlProp2DAcoTTIDenQ_DEO2_FDTD(; b = b, θ = θ, ntrec = ntrec, dz = dz, dx = dx, 
    dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,1] .= 1500;
m₀[:,:,2] .= 0.1;
m₀[:,:,3] .= 0.2;
J = jacobian(F, m₀)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Construct and apply the linearized Jacobian operator (method 2)

For this example we assume in the model that v, ϵ, η are active parameters, and b, θ are passive.

  1. create the model vector and set parameters v, ϵ, η
  2. create the Jacobian operator J at point m₀ directly by construction of a jet.
  3. create a random model perturbation vector δm.
  4. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  5. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm. 
m₀ = zeros(Float32, nz, nx, 3);
m₀[:,:,1] .= 1500;
m₀[:,:,2] .= 0.1;
m₀[:,:,3] .= 0.2;
J = JopLnProp2DAcoTTIDenQ_DEO2_FDTD(; m₀ = m₀, b = b, θ = θ, ntrec = ntrec, dz = dz, dx = dx, 
    dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Required Parameters

  • m₀ the point at which the jet is linearized. Note this argument is required in the constuctor for JopLnProp2DAcoTTIDenQ_DEO2_FDTD but not JopNlProp2DAcoTTIDenQ_DEO2_FDTD. This constuctor is shown in the last example above. Please note that you must consider which parameters are active and passive, per the discussion on model parameters above and examples.
  • dtmod The sample rate for the modeled data. You can establish a lower bound for the modeling sample rate with the following expression: dt = 0.75 * 0.38 * max(dx,dz) / maximum(v). Note that the usual Courant–Friedrichs–Lewy condition (CFL condition) for stability in finite difference modeling is modified heuristically to include the impact of the visco-acoustic implementation of this operator, requiring a 25% smaller smaple rate.
  • dtrec The number of time samples in the recorded data. Note requirement that dtrec be an even multiple of dtmod.
  • ntrec The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined from ntrec, dtrec, and dtmod.

Optional Parameters

Defaults for arguments are shown inside square brackets.

  • b [ones(Float32,0,0)] the buoyancy (reciprocal density) array.
  • θ [zeros(Float32,0,0)] the symmetry axies tilt angle in radians. Assumed to be zero if not specified. θ must be an array the same size as buoyancy b.
  • f [0.85] See the discussion above regarding the Pseudo-acoustic approximation used here
  • srcfieldfile [joinpath(tempdir(), "field-55e40154-73c6-434c-b029-1fe5243e0e7d.bin")] the full path to a scratch file used for the serializationof the compressed nonlinear source wavefield.
  • comptype [nothing] the type of compression to use for the serialization of the nonlinear source wavefield. The type of compression must be one of:
    • nothing - no compression.
    • Float32 - if wavefield is Float64, then do a simple conversion to Float32.
    • UInt32 - compression using CvxCompress (windowing + 2D wavelet transform + thresholding + quantization + run-length-encoding).
  • compscale [1e-2] determines the thresholding for the compression of the nonlinear source wavefield prior to serialization. Larger values mean more aggressive compression. You can likely increase compscale to 1.0 before you start to notice differences in output from Jacobian operations.
  • nz_subcube, nx_subcube [32] The Z and X sizes of windows used for compression of the nonlinear source wavefield with CvxCompress. Note the requirement [8 <= n*_subcube <=256].
  • isinterior [false] boolean flag that indicates how the nonlinear source wavefield is saved. For large models, operation will be faster with isinterior = true, but the linearization correctness test may fail.
    • true the entire model including absorbing boundaries is serialized and deserialized
    • false the interior part of the model excluding absorbing boundaries is serialized and deserialized
  • sz [0.0] Array of source Z coordinate. Note that if multiple sources are provided, they will will be injected simultaneously during finite difference evolution.
  • sx [0.0] Array of source X coordinate.
  • st [0.0] Array of source delay times.
  • interpmethod [:hicks] Type of physical interpolation for sources and receivers. For locations that are not on the physical grid coordinates, interpolation is used either to inject or extract information. interpmethod must be one of:
    • :hicks Hicks 2D sinc interpolation (up to 8x8 nonzero points per location)
    • :linear bilinear interpolation (up to 2x2 nonzero points per location)
  • rz [[0.0]] 2D array of receiver Z coordinates
  • rx [[0.0]] 2D array of receiver Z coordinates
  • z0 [0.0] Origin of physical coordinates for the Z dimension.
  • x0 [0.0] Origin of physical coordinates for the X dimension.
  • dz [10.0] Spacing of physical coordinates in the Z dimension.
  • dx [10.0] Spacing of physical coordinates in the X dimension.
  • freqQ [5.0] The center frequency for the Maxwell body approximation to dissipation only attenuation. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qMin [0.1] The minimum value for Qp at the boundary of the model used in our Maxwell body approximation to dissipation only attenuation. This is not a physically meaningful value for Qp, as we use the attenuation to implement absorbing boundary conditions and eliminate outgoing waves on the boundaries of the computational domain. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qInterior [100.0] the value for Qp in the interior of the model used in our Maxwell body approximation to dissipation only attenuation. This is the value for Qp away from the absorbing boundaries and is a physically meaningful value. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • padz, padx [0.0], [0.0] - apply extra padding to the survey determined aperture in Ginsu. Please see Ginsu for more information.
  • nbz_cache, nbx_cache [512], [8] The size of cache blocks in the Z and X dimensions. In general the cache block in the Z (fast) dimension should be ≥ the entire size of that dimension, and the cache block size in the slower dimension is generally small in order to allow the entire block to fit in cache.
  • nsponge [50] The number of grid cells to use for the absorbing boundary. For high fidelity modeling this should be > 60 grid points, but can be significantly smaller for some use cases like low frequency full waveform inversion.
  • wavelet [WaveletCausalRicker(f=5.0)] The source wavelet, can be specified as either a Wavelet type or an array.
  • freesurface [false] Determines if a free surface (true) or absorbing (false) top boundary condition is applied.
  • imgcondition ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint only exists for "standard" imaging condition currently.
  • RTM_weight [0.5] determines the balance of short wavelengths and long wavelengths in the imaging condition. A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition.
  • nthreads [Sys.CPU_THREADS] The number of threads to use for OpenMP parallelization of the modeling.
  • reportinterval [500] The interval at which information about the propagtion is logged.

See also: Ginsu, WaveletSine, WaveletRicker, WaveletMinPhaseRicker, WaveletDerivRicker, WaveletCausalRicker, WaveletOrmsby, WaveletMinPhaseOrmsby

source
JetPackWaveFD.JopNlProp2DAcoVTIDenQ_DEO2_FDTD — Method
JopNlProp2DAcoVTIDenQ_DEO2_FDTD(; kwargs...)
JopLnProp2DAcoVTIDenQ_DEO2_FDTD(; m₀, kwargs...)

Create a Jets nonlinear or linearized operator for 2D pseudo- visco-acoustic, vertical transverse isotropy, variable density modeling.

Model Parameters

This propagator operates with four model parameters, as shown in the table below.

ParameterDescription
vP wave velocity
ϵModified Thomsen's weak anisotropy parameter
ηModified Alkhalifah's weak anisotropy parameter η
bbuoyancy (reciprocal density)

Pseudo-acoustic approximation

With the pseudo-acoustic approximation we employ, the transformation to self-adjoint form requires a parameter f representing the average ratio of shear-wave to p-wave velocity, and a modfication to Alkhalifah's weak anisotropy parameter η. We show formulas for these two parameters below. The parameter f is specified as a scalar for the entire model, with default value 0.85, implying a shear velocity around 38% of the p wave velocity.

  • η = sqrt[2(ϵ-δ)/(f+2ϵ)]
  • f = 1 - Vₛ² / Vₚ²

For more information about the pseudo-acoustic approximation we employ please see https://library.seg.org/doi/10.1190/segam2016-13878451.1.

Active and passive model parameters

There are two modes of operation that define different sets of active and passive parameters, velocity only and velocity and anisotropy.

An active parameter can be inverted for using the Jacobian linearized operator machinery, and a passive parameter is constant. The two modes are:

ModeActive ParametersPassive Parameters
velocity only[v][ϵ, η, b]
velocity and anisotropy[v, ϵ, η][b]

To make a parameter passive, you pass the array for that parameter in the constructor for the operator, implying that it is in state and does not change with the action of the operator. Parameters that are not passed to the constructor are assumed to be active, and will be part of the model that the operator acts on, stored as a 3D array with the following indices: [Z coord][X coord][Active Parameter]

Examples

Model and acquisition geometry setup

  1. load modules Jets, WaveFD, and JetPackWaveFD
  2. set up the model discretization, coordinate size and spacing
  3. set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
  4. create constant models for v, ϵ, η, b
using Jets, WaveFD, JetPackWaveFD
nz,nx = 200, 100                      # spatial discretization size
dz,dx = 20.0, 20.0                    # spatial discretization sampling
ntrec = 1101                          # number of temporal samples in recorded data
dtrec = 0.0100                        # sample rate for recorded data
dtmod = 0.0025                        # sample rate for modeled data
wavelet = WaveletCausalRicker(f=5.0)  # type of wavelet to use (source signature) 
sz = dz                               # source Z location 
sx = dx*(nx/2)                        # source X location 
rx = dx*[0:0.5:nx-1;];                # array of receiver X locations
rz = 2*dz*ones(length(0:0.5:nx-1));   # array of receiver Z locations
b = ones(Float32, nz, nx);            # constant buoyancy
ϵ = 0.1 * ones(Float32, nz, nx);      # constant epsilon
η = 0.2 * ones(Float32, nz, nx);      # constant eta

Construct and apply the nonlinear operator (v is active; ϵ, η, b are passive)

  1. create the nonlinear operator F
  2. create the model vector and set parameter v
  3. perform nonlinear forward modeling with model m₀ and return the resulting modeled data in d
F = JopNlProp2DAcoVTIDenQ_DEO2_FDTD(; b = b, ϵ = ϵ, η = η, ntrec = ntrec, dtrec = dtrec, 
    dtmod = dtmod, dz = dz, dx = dx, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,1] .= 1500;
d = F*m₀;              # forward nonlinear op

Construct and apply the nonlinear operator (v, ϵ, η are active; b is passive)

  1. create the nonlinear operator F
  2. create the model vector and set parameters v, ϵ, η
  3. perform nonlinear forward modeling with model m₀ and return the resulting modeled data in d
F = JopNlProp2DAcoVTIDenQ_DEO2_FDTD(; b = b, ntrec = ntrec, dtrec = dtrec, dtmod = dtmod,
    dz = dz, dx = dx, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,1] .= 1500;
m₀[:,:,2] .= 0.1;
m₀[:,:,3] .= 0.2;
d = F*m₀;              # forward nonlinear op

Construct and apply the linearized Jacobian operator (method 1)

For this example we assume in the model that v, ϵ, η are active parameters, and b is passive.

  1. create the nonlinear operator F directly by construction of a jet.
  2. create the model vector and set parameters v, ϵ, η
  3. create the Jacobian operator J by linearization of F at point m₀
  4. create a random model perturbation vector δm.
  5. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  6. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm.

Note that the Jacobian operators J and J' require the serialized nonlinear forward wavefield, and this is generated automatically whenever required. If you watch the logging to standard out from this example, you will first see the finite difference evolution for the nonlinear forward, followed by the linearized forward, and finally the linearized adjoint.

F = JopNlProp2DAcoVTIDenQ_DEO2_FDTD(; b = b, ntrec = ntrec, dz = dz, dx = dx, 
    dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,1] .= 1500;
m₀[:,:,2] .= 0.1;
m₀[:,:,3] .= 0.2;
J = jacobian(F, m₀)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Construct and apply the linearized Jacobian operator (method 2)

For this example we assume in the model that v, ϵ, η are active parameters, and b is passive.

  1. create the model vector and set parameters v, ϵ, η
  2. create the Jacobian operator J at point m₀ directly by construction of a jet.
  3. create a random model perturbation vector δm.
  4. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  5. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm. 
m₀ = zeros(Float32, nz, nx, 3);
m₀[:,:,1] .= 1500;
m₀[:,:,2] .= 0.1;
m₀[:,:,3] .= 0.2;
J = JopLnProp2DAcoVTIDenQ_DEO2_FDTD(; m₀ = m₀, b = b, ntrec = ntrec, dz = dz, dx = dx, 
    dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Required Parameters

  • m₀ the point at which the jet is linearized. Note this argument is required in the constuctor for JopLnProp2DAcoVTIDenQ_DEO2_FDTD but not JopNlProp2DAcoVTIDenQ_DEO2_FDTD. This constuctor is shown in the last example above. Please note that you must consider which parameters are active and passive, per the discussion on model parameters above and examples.
  • dtmod The sample rate for the modeled data. You can establish a lower bound for the modeling sample rate with the following expression: dt = 0.75 * 0.38 * max(dx,dz) / maximum(v). Note that the usual Courant–Friedrichs–Lewy condition (CFL condition) for stability in finite difference modeling is modified heuristically to include the impact of the visco-acoustic implementation of this operator, requiring a 25% smaller smaple rate.
  • dtrec The number of time samples in the recorded data. Note requirement that dtrec be an even multiple of dtmod.
  • ntrec The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined from ntrec, dtrec, and dtmod.

Optional Parameters

Defaults for arguments are shown inside square brackets.

  • b [ones(Float32,0,0)] the buoyancy (reciprocal density) array.
  • f [0.85] See the discussion above regarding the Pseudo-acoustic approximation used here
  • srcfieldfile [joinpath(tempdir(), "field-24b82d11-9647-4676-a2e8-c34c19b03b99.bin")] the full path to a scratch file used for the serializationof the compressed nonlinear source wavefield.
  • comptype [nothing] the type of compression to use for the serialization of the nonlinear source wavefield. The type of compression must be one of:
    • nothing - no compression.
    • Float32 - if wavefield is Float64, then do a simple conversion to Float32.
    • UInt32 - compression using CvxCompress (windowing + 2D wavelet transform + thresholding + quantization + run-length-encoding).
  • compscale [1e-2] determines the thresholding for the compression of the nonlinear source wavefield prior to serialization. Larger values mean more aggressive compression. You can likely increase compscale to 1.0 before you start to notice differences in output from Jacobian operations.
  • nz_subcube, nx_subcube [32] The Z and X sizes of windows used for compression of the nonlinear source wavefield with CvxCompress. Note the requirement [8 <= n*_subcube <=256].
  • isinterior [false] boolean flag that indicates how the nonlinear source wavefield is saved. For large models, operation will be faster with isinterior = true, but the linearization correctness test may fail.
    • true the entire model including absorbing boundaries is serialized and deserialized
    • false the interior part of the model excluding absorbing boundaries is serialized and deserialized
  • sz [0.0] Array of source Z coordinate. Note that if multiple sources are provided, they will will be injected simultaneously during finite difference evolution.
  • sx [0.0] Array of source X coordinate.
  • st [0.0] Array of source delay times.
  • interpmethod [:hicks] Type of physical interpolation for sources and receivers. For locations that are not on the physical grid coordinates, interpolation is used either to inject or extract information. interpmethod must be one of:
    • :hicks Hicks 2D sinc interpolation (up to 8x8 nonzero points per location)
    • :linear bilinear interpolation (up to 2x2 nonzero points per location)
  • rz [[0.0]] 2D array of receiver Z coordinates
  • rx [[0.0]] 2D array of receiver Z coordinates
  • z0 [0.0] Origin of physical coordinates for the Z dimension.
  • x0 [0.0] Origin of physical coordinates for the X dimension.
  • dz [10.0] Spacing of physical coordinates in the Z dimension.
  • dx [10.0] Spacing of physical coordinates in the X dimension.
  • freqQ [5.0] The center frequency for the Maxwell body approximation to dissipation only attenuation. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qMin [0.1] The minimum value for Qp at the boundary of the model used in our Maxwell body approximation to dissipation only attenuation. This is not a physically meaningful value for Qp, as we use the attenuation to implement absorbing boundary conditions and eliminate outgoing waves on the boundaries of the computational domain. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qInterior [100.0] the value for Qp in the interior of the model used in our Maxwell body approximation to dissipation only attenuation. This is the value for Qp away from the absorbing boundaries and is a physically meaningful value. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • padz, padx [0.0], [0.0] - apply extra padding to the survey determined aperture in Ginsu. Please see Ginsu for more information.
  • nbz_cache, nbx_cache [512], [8] The size of cache blocks in the Z and X dimensions. In general the cache block in the Z (fast) dimension should be ≥ the entire size of that dimension, and the cache block size in the slower dimension is generally small in order to allow the entire block to fit in cache.
  • nsponge [50] The number of grid cells to use for the absorbing boundary. For high fidelity modeling this should be > 60 grid points, but can be significantly smaller for some use cases like low frequency full waveform inversion.
  • wavelet [WaveletCausalRicker(f=5.0)] The source wavelet, can be specified as either a Wavelet type or an array.
  • freesurface [false] Determines if a free surface (true) or absorbing (false) top boundary condition is applied.
  • imgcondition ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint only exists for "standard" imaging condition currently.
  • RTM_weight [0.5] determines the balance of short wavelengths and long wavelengths in the imaging condition. A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition.
  • nthreads [Sys.CPU_THREADS] The number of threads to use for OpenMP parallelization of the modeling.
  • reportinterval [500] The interval at which information about the propagtion is logged.

See also: Ginsu, WaveletSine, WaveletRicker, WaveletMinPhaseRicker, WaveletDerivRicker, WaveletCausalRicker, WaveletOrmsby, WaveletMinPhaseOrmsby

source
JetPackWaveFD.JopNlProp3DAcoIsoDenQ_DEO2_FDTD — Method
JopNlProp3DAcoIsoDenQ_DEO2_FDTD(; kwargs...)
JopLnProp3DAcoIsoDenQ_DEO2_FDTD(; v, kwargs...)

Create a Jets nonlinear or linearized operator for 3D visco-acoustic, isotropic, variable density modeling.

Model Parameters

This propagator operates with two model parameters, as shown in the table below.

ParameterDescription
vP wave velocity
bbuoyancy (reciprocal density)

Velocity v and buoyancy b can be active parameters and inverted for using the Jacobian machinery, or passive parameters that are constant. If a parameter is active there will be a component of model assigned to it.

The model is 4D with the slowest dimension for model parameter. The earth model property maps to the slow- fimension array index via the state(F).active_modelset dictionary.

Examples

Model and acquisition geometry setup

  1. load modules Jets, WaveFD, and JetPackWaveFD
  2. set up the model discretization, coordinate size and spacing
  3. set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
  4. create constant buoyancy array
using Jets, WaveFD, JetPackWaveFD
nz,ny,nx = 100, 80, 60                        # spatial discretization size
dz,dy,dx = 20.0, 20.0, 20.0                   # spatial discretization sampling
ntrec = 1101                                  # number of temporal samples in recorded data
dtrec = 0.0100                                # sample rate for recorded data
dtmod = 0.0025                                # sample rate for modeled data
wavelet = WaveletCausalRicker(f=5.0)          # type of wavelet to use (source signature) 
sz = dz                                       # source Z location 
sy = dy*(ny/2)                                # source Y location 
sx = dx*(nx/2)                                # source X location 
rz = [dz for iy = 1:ny, ix = 1:nx][:];        # Array of receiver Y locations 
ry = [(iy-1)*dy for iy = 1:ny, ix = 1:nx][:]; # Array of receiver Y locations
rx = [(ix-1)*dx for iy = 1:ny, ix=1:nx][:];   # Array of receiver X locations
b = ones(Float32, nz, ny, nx);                # buoyancy model (reciprocal density)

Construct and apply the nonlinear operator

  1. create the nonlinear operator F
  2. create the constant velocity model m₀
  3. perform nonlinear forward modeling with constant velocity model m₀ and return the resulting modeled data in d
F = JopNlProp3DAcoIsoDenQ_DEO2_FDTD(; b = b, isinterior=true, nsponge = 10, ntrec = ntrec, 
    dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = 1500 .* ones(domain(F));
d = F*m₀;              # forward nonlinear op

Construct and apply the linearized Jacobian operator (method 1)

  1. create the nonlinear operator F directly by construction of a jet.
  2. create the constant velocity model m₀
  3. create the Jacobian operator J by linearization of F at point m₀
  4. create a random model perturbation vector δm.
  5. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  6. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm.

Note that the Jacobian operators J and J' require the serialized nonlinear forward wavefield, and this is generated automatically whenever required. If you watch the logging to standard out from this example, you will first see the finite difference evolution for the nonlinear forward, followed by the linearized forward, and finally the linearized adjoint.

F = JopNlProp3DAcoIsoDenQ_DEO2_FDTD(; b = b, isinterior=true, nsponge = 10, ntrec = ntrec, 
    dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = 1500 .* ones(domain(F));
J = jacobian(F, m₀)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Construct and apply the linearized Jacobian operator (method 2)

  1. create the constant velocity model m₀
  2. create the Jacobian operator J at point m₀ directly by construction of a jet.
  3. create a random model perturbation vector δm.
  4. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  5. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm. 
m₀ = 1500 .* ones(Float32, nz, ny, nx);       # velocity model
J = JopLnProp3DAcoIsoDenQ_DEO2_FDTD(; v = m₀, b = b, isinterior=true, nsponge = 10, ntrec = ntrec, 
    dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sy = sy, sx = sx, ry = ry, rz = rz, rx = rx)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Required Parameters

  • v the point at which the jet is linearized. Note this argument is required in the constuctor for JopLnProp3DAcoIsoDenQ_DEO2_FDTD but not JopNlProp3DAcoIsoDenQ_DEO2_FDTD. This constuctor is shown in the last example above.
  • dtmod The sample rate for the modeled data. You can establish a lower bound for the modeling sample rate with the following expression: dt = 0.75 * 0.38 * max(dx,dz) / maximum(v). Note that the usual Courant–Friedrichs–Lewy condition (CFL condition) for stability in finite difference modeling is modified heuristically to include the impact of the visco-acoustic implementation of this operator, requiring a 25% smaller smaple rate.
  • dtrec The number of time samples in the recorded data. Note requirement that dtrec be an even multiple of dtmod.
  • ntrec The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined from ntrec, dtrec, and dtmod.

Optional Parameters

Defaults for arguments are shown inside square brackets.

  • b [ones(Float32,0,0,0)] the buoyancy (reciprocal density) array.
  • srcfieldfile [joinpath(tempdir(), "field-977465bb-cd8e-46ef-91e2-49b567797c0b.bin")] the full path to a scratch file used for the serializationof the compressed nonlinear source wavefield.
  • comptype [nothing] the type of compression to use for the serialization of the nonlinear source wavefield. The type of compression must be one of:
    • nothing - no compression.
    • Float32 - if wavefield is Float64, then do a simple conversion to Float32.
    • UInt32 - compression using CvxCompress (windowing + 2D wavelet transform + thresholding + quantization + run-length-encoding).
  • compscale [1e-2] determines the thresholding for the compression of the nonlinear source wavefield prior to serialization. Larger values mean more aggressive compression. You can likely increase compscale to 1.0 before you start to notice differences in output from Jacobian operations.
  • nz_subcube, ny_subcube, nx_subcube [32] The Z, Y, and X sizes of windows used for compression of the nonlinear source wavefield with CvxCompress. Note the requirement [8 <= n*_subcube <=256].
  • isinterior [false] boolean flag that indicates how the nonlinear source wavefield is saved. For large models, operation will be faster with isinterior = true, but the linearization correctness test may fail.
    • true the entire model including absorbing boundaries is serialized and deserialized
    • false the interior part of the model excluding absorbing boundaries is serialized and deserialized
  • sz [0.0] Array of source Z coordinate. Note that if multiple sources are provided, they will will be injected simultaneously during finite difference evolution.
  • sy [0.0] Array of source Y coordinate.
  • sx [0.0] Array of source X coordinate.
  • st [0.0] Array of source delay times.
  • interpmethod [:hicks] Type of physical interpolation for sources and receivers. For locations that are not on the physical grid coordinates, interpolation is used either to inject or extract information. interpmethod must be one of:
    • :hicks Hicks 3D sinc interpolation (up to 8x8x8 nonzero points per location)
    • :linear bilinear interpolation (up to 2x2x2 nonzero points per location)
  • rz [[0.0]] 2D array of receiver Z coordinates
  • ry [[0.0]] 2D array of receiver Y coordinates
  • rx [[0.0]] 2D array of receiver Z coordinates
  • z0 [0.0] Origin of physical coordinates for the Z dimension.
  • y0 [0.0] Origin of physical coordinates for the Y dimension.
  • x0 [0.0] Origin of physical coordinates for the X dimension.
  • dz [10.0] Spacing of physical coordinates in the Z dimension.
  • dy [10.0] Spacing of physical coordinates in the Y dimension.
  • dx [10.0] Spacing of physical coordinates in the X dimension.
  • freqQ [5.0] The center frequency for the Maxwell body approximation to dissipation only attenuation. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qMin [0.1] The minimum value for Qp at the boundary of the model used in our Maxwell body approximation to dissipation only attenuation. This is not a physically meaningful value for Qp, as we use the attenuation to implement absorbing boundary conditions and eliminate outgoing waves on the boundaries of the computational domain. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qInterior [100.0] the value for Qp in the interior of the model used in our Maxwell body approximation to dissipation only attenuation. This is the value for Qp away from the absorbing boundaries and is a physically meaningful value. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • padz, pady, padx [0.0], [0.0], [0.0] - apply extra padding to the survey determined aperture in Ginsu. Please see Ginsu for more information.
  • nbz_cache, nby_cache, nbx_cache [512], [8], [8] The size of cache blocks in the Z, X, and Y dimensions. In general the cache block in the Z (fast) dimension should be ≥ the entire size of that dimension, and the cache block size in the slower dimensions is generally small in order to allow the entire block to fit in cache.
  • nsponge [50] The number of grid cells to use for the absorbing boundary. For high fidelity modeling this should be > 60 grid points, but can be significantly smaller for some use cases like low frequency full waveform inversion.
  • wavelet [WaveletCausalRicker(f=5.0)] The source wavelet, can be specified as either a Wavelet type or an array.
  • freesurface [false] Determines if a free surface (true) or absorbing (false) top boundary condition is applied.
  • imgcondition ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint only exists for "standard" imaging condition currently.
  • RTM_weight [0.5] determines the balance of short wavelengths and long wavelengths in the imaging condition. A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition.
  • nthreads [Sys.CPU_THREADS] The number of threads to use for OpenMP parallelization of the modeling.
  • reportinterval [500] The interval at which information about the propagtion is logged.

See also: Ginsu, WaveletSine, WaveletRicker, WaveletMinPhaseRicker, WaveletDerivRicker, WaveletCausalRicker, WaveletOrmsby, WaveletMinPhaseOrmsby

source
JetPackWaveFD.JopNlProp3DAcoTTIDenQ_DEO2_FDTD — Method
JopNlProp3DAcoTTIDenQ_DEO2_FDTD(; kwargs...)
JopLnProp3DAcoTTIDenQ_DEO2_FDTD(; m₀, kwargs...)

Create a Jets nonlinear or linearized operator for 3D pseudo- visco-acoustic, tilted transverse isotropy, variable density modeling.

Model Parameters

This propagator operates with four model parameters, as shown in the table below.

ParameterDescription
vP wave velocity
ϵModified Thomsen's weak anisotropy parameter
ηModified Alkhalifah's weak anisotropy parameter η
bbuoyancy (reciprocal density)
θsymmetry axis tilt angle from vertical (radians)
ϕsymmetry axis aziumuth angle CCW from x (radians)

Pseudo-acoustic approximation

With the pseudo-acoustic approximation we employ, the transformation to self-adjoint form requires a parameter f representing the average ratio of shear-wave to p-wave velocity, and a modfication to Alkhalifah's weak anisotropy parameter η. We show formulas for these two parameters below. The parameter f is specified as a scalar for the entire model, with default value 0.85, implying a shear velocity around 38% of the p wave velocity.

  • η = sqrt[2(ϵ-δ)/(f+2ϵ)]
  • f = 1 - Vₛ² / Vₚ²

For more information about the pseudo-acoustic approximation we employ please see https://library.seg.org/doi/10.1190/segam2016-13878451.1.

Active and passive model parameters

There are two modes of operation that define different sets of active and passive parameters, velocity only and velocity and anisotropy.

An active parameter can be inverted for using the Jacobian linearized operator machinery, and a passive parameter is constant. The two modes are:

ModeActive ParametersPassive Parameters
velocity only[v][ϵ, η, b, θ, ϕ ]
velocity and anisotropy[v, ϵ, η][b, θ, ϕ ]

To make a parameter passive, you pass the array for that parameter in the constructor for the operator, implying that it is in state and does not change with the action of the operator. Parameters that are not passed to the constructor are assumed to be active, and will be part of the model that the operator acts on, stored as a 3D array with the following indices: [Z coord][X coord][Active Parameter]

Examples

Model and acquisition geometry setup

  1. load modules Jets, WaveFD, and JetPackWaveFDFDFD
  2. set up the model discretization, coordinate size and spacing
  3. set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
  4. create constant models for v, ϵ, η, b, θ, ϕ
using Jets, WaveFD, JetPackWaveFDFDFD
nz,ny,nx = 100, 80, 60                        # spatial discretization size
dz,dy,dx = 20.0, 20.0, 20.0                   # spatial discretization sampling
ntrec = 1101                                  # number of temporal samples in recorded data
dtrec = 0.0100                                # sample rate for recorded data
dtmod = 0.0025                                # sample rate for modeled data
wavelet = WaveletCausalRicker(f=5.0)          # type of wavelet to use (source signature) 
sz = dz                                       # source Z location 
sy = dy*(ny/2)                                # source Y location 
sx = dx*(nx/2)                                # source X location 
rz = [dz for iy = 1:ny, ix = 1:nx][:];        # Array of receiver Y locations 
ry = [(iy-1)*dy for iy = 1:ny, ix = 1:nx][:]; # Array of receiver Y locations
rx = [(ix-1)*dx for iy = 1:ny, ix=1:nx][:];   # Array of receiver X locations
b = ones(Float32, nz, ny, nx);                # constant buoyancy
ϵ = 0.1 * ones(Float32, nz, ny, nx);          # constant epsilon
η = 0.2 * ones(Float32, nz, ny, nx);          # constant eta
θ = (π/8) * ones(Float32, nz, ny, nx);        # constant tilt angle
ϕ = (π/3) * ones(Float32, nz, ny, nx);        # constant azimuth angle

Construct and apply the nonlinear operator (v is active; ϵ, η, b, θ, ϕ are passive)

  1. create the nonlinear operator F
  2. create the constant velocity model m₀
  3. perform nonlinear forward modeling with constant velocity model m₀ and return the resulting modeled data in d
F = JopNlProp3DAcoTTIDenQ_DEO2_FDTD(; b = b, ϵ = ϵ, η = η, θ = θ, ϕ = ϕ, isinterior=true, nsponge = 10, 
    ntrec = ntrec, dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, 
    wavelet = wavelet, sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,:,1] .= 1500;
d = F*m₀;              # forward nonlinear op

Construct and apply the nonlinear operator (v, ϵ, η are active; b, θ, ϕ are passive)

  1. create the nonlinear operator F
  2. create the model vector and set parameters v, ϵ, η
  3. perform nonlinear forward modeling with model m₀ and return the resulting modeled data in d
F = JopNlProp3DAcoTTIDenQ_DEO2_FDTD(; b = b, θ = θ, ϕ = ϕ, isinterior=true, nsponge = 10, 
    ntrec = ntrec, dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, 
    wavelet = wavelet, sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,:,1] .= 1500;
m₀[:,:,:,2] .= 0.1;
m₀[:,:,:,3] .= 0.2;
d = F*m₀;              # forward nonlinear op

Construct and apply the linearized Jacobian operator (method 1)

For this example we assume in the model that v, ϵ, η are active parameters, and b, θ, ϕ are passive.

  1. create the nonlinear operator F directly by construction of a jet.
  2. create the model vector and set parameters v, ϵ, η
  3. create the Jacobian operator J by linearization of F at point m₀
  4. create a random model perturbation vector δm.
  5. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  6. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm.

Note that the Jacobian operators J and J' require the serialized nonlinear forward wavefield, and this is generated automatically whenever required. If you watch the logging to standard out from this example, you will first see the finite difference evolution for the nonlinear forward, followed by the linearized forward, and finally the linearized adjoint.

F = JopNlProp3DAcoTTIDenQ_DEO2_FDTD(; b = b, θ = θ, ϕ = ϕ, isinterior=true, nsponge = 10, 
    ntrec = ntrec, dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,:,1] .= 1500;
m₀[:,:,:,2] .= 0.1;
m₀[:,:,:,3] .= 0.2;
J = jacobian(F, m₀)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Construct and apply the linearized Jacobian operator (method 2)

For this example we assume in the model that v, ϵ, η are active parameters, and b, θ, ϕ are passive.

  1. create the constant velocity model m₀
  2. create the Jacobian operator J at point m₀ directly by construction of a jet.
  3. create a random model perturbation vector δm.
  4. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  5. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm. 
m₀ = zeros(Float32, nz, ny, nx, 3);
m₀[:,:,:,1] .= 1500;
m₀[:,:,:,2] .= 0.1;
m₀[:,:,:,3] .= 0.2;
J = JopLnProp3DAcoTTIDenQ_DEO2_FDTD(; m₀ = m₀, b = b, θ = θ, ϕ = ϕ, isinterior=true, nsponge = 10, 
    ntrec = ntrec, dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, 
    wavelet = wavelet, sz = sz, sy = sy, sx = sx, ry = ry, rz = rz, rx = rx)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Required Parameters

  • m₀ the point at which the jet is linearized. Note this argument is required in the constuctor for JopLnProp2DAcoTTIDenQ_DEO2_FDTD but not JopNlProp2DAcoTTIDenQ_DEO2_FDTD. This constuctor is shown in the last example above. Please note that you must consider which parameters are active and passive, per the discussion on model parameters above and examples.
  • dtmod The sample rate for the modeled data. You can establish a lower bound for the modeling sample rate with the following expression: dt = 0.75 * 0.38 * max(dx,dz) / maximum(v). Note that the usual Courant–Friedrichs–Lewy condition (CFL condition) for stability in finite difference modeling is modified heuristically to include the impact of the visco-acoustic implementation of this operator, requiring a 25% smaller smaple rate.
  • dtrec The number of time samples in the recorded data. Note requirement that dtrec be an even multiple of dtmod.
  • ntrec The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined from ntrec, dtrec, and dtmod.

Optional Parameters

Defaults for arguments are shown inside square brackets.

  • b [ones(Float32,0,0,0)] the buoyancy (reciprocal density) array.
  • ϕ [zeros(Float32,0,0,0)] the symmetry axies azimuth angle counter-clockwise from the x axis, in radians. Assumed to be zero if not specified.
  • θ [zeros(Float32,0,0,0)] the symmetry axies tilt angle in radians. Assumed to be zero if not specified. θ must be an array the same size as buoyancy b.
  • f [0.85] See the discussion above regarding the Pseudo-acoustic approximation used here
  • srcfieldfile ["field-4fee5977-1631-499f-9f65-93b02219bb06.bin"] the full path to a scratch file used for the serializationof the compressed nonlinear source wavefield.
  • comptype [nothing] the type of compression to use for the serialization of the nonlinear source wavefield. The type of compression must be one of:
    • nothing - no compression.
    • Float32 - if wavefield is Float64, then do a simple conversion to Float32.
    • UInt32 - compression using CvxCompress (windowing + 2D wavelet transform + thresholding + quantization + run-length-encoding).
  • compscale [1e-2] determines the thresholding for the compression of the nonlinear source wavefield prior to serialization. Larger values mean more aggressive compression. You can likely increase compscale to 1.0 before you start to notice differences in output from Jacobian operations.
  • nz_subcube, ny_subcube, nx_subcube [32] The Z, Y, and X sizes of windows used for compression of the nonlinear source wavefield with CvxCompress. Note the requirement [8 <= n*_subcube <=256].
  • isinterior [false] boolean flag that indicates how the nonlinear source wavefield is saved. For large models, operation will be faster with isinterior = true, but the linearization correctness test may fail.
    • true the entire model including absorbing boundaries is serialized and deserialized
    • false the interior part of the model excluding absorbing boundaries is serialized and deserialized
  • sz [0.0] Array of source Z coordinate. Note that if multiple sources are provided, they will will be injected simultaneously during finite difference evolution.
  • sy [0.0] Array of source Y coordinate.
  • sx [0.0] Array of source X coordinate.
  • st [0.0] Array of source delay times.
  • interpmethod [:hicks] Type of physical interpolation for sources and receivers. For locations that are not on the physical grid coordinates, interpolation is used either to inject or extract information. interpmethod must be one of:
    • :hicks Hicks 3D sinc interpolation (up to 8x8x8 nonzero points per location)
    • :linear bilinear interpolation (up to 2x2x2 nonzero points per location)
  • rz [[0.0]] 2D array of receiver Z coordinates
  • ry [[0.0]] 2D array of receiver Y coordinates
  • rx [[0.0]] 2D array of receiver Z coordinates
  • z0 [0.0] Origin of physical coordinates for the Z dimension.
  • y0 [0.0] Origin of physical coordinates for the Y dimension.
  • x0 [0.0] Origin of physical coordinates for the X dimension.
  • dz [10.0] Spacing of physical coordinates in the Z dimension.
  • dy [10.0] Spacing of physical coordinates in the Y dimension.
  • dx [10.0] Spacing of physical coordinates in the X dimension.
  • freqQ [5.0] The center frequency for the Maxwell body approximation to dissipation only attenuation. Please see JetPackWaveFDFDFD package documentation for more information concerning the attenuation model.
  • qMin [0.1] The minimum value for Qp at the boundary of the model used in our Maxwell body approximation to dissipation only attenuation. This is not a physically meaningful value for Qp, as we use the attenuation to implement absorbing boundary conditions and eliminate outgoing waves on the boundaries of the computational domain. Please see JetPackWaveFDFDFD package documentation for more information concerning the attenuation model.
  • qInterior [100.0] the value for Qp in the interior of the model used in our Maxwell body approximation to dissipation only attenuation. This is the value for Qp away from the absorbing boundaries and is a physically meaningful value. Please see JetPackWaveFDFDFD package documentation for more information concerning the attenuation model.
  • padz, pady, padx [0.0], [0.0], [0.0] - apply extra padding to the survey determined aperture in Ginsu. Please see Ginsu for more information.
  • nbz_cache, nby_cache, nbx_cache [512], [8], [8] The size of cache blocks in the Z, X, and Y dimensions. In general the cache block in the Z (fast) dimension should be ≥ the entire size of that dimension, and the cache block size in the slower dimensions is generally small in order to allow the entire block to fit in cache.
  • nsponge [50] The number of grid cells to use for the absorbing boundary. For high fidelity modeling this should be > 60 grid points, but can be significantly smaller for some use cases like low frequency full waveform inversion.
  • wavelet [WaveletCausalRicker(f=5.0)] The source wavelet, can be specified as either a Wavelet type or an array.
  • freesurface [false] Determines if a free surface (true) or absorbing (false) top boundary condition is applied.
  • imgcondition ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint only exists for "standard" imaging condition currently.
  • RTM_weight [0.5] determines the balance of short wavelengths and long wavelengths in the imaging condition. A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition.
  • nthreads [Sys.CPU_THREADS] The number of threads to use for OpenMP parallelization of the modeling.
  • reportinterval [500] The interval at which information about the propagtion is logged.

See also: Ginsu, WaveletSine, WaveletRicker, WaveletMinPhaseRicker, WaveletDerivRicker, WaveletCausalRicker, WaveletOrmsby, WaveletMinPhaseOrmsby

source
JetPackWaveFD.JopNlProp3DAcoVTIDenQ_DEO2_FDTD — Method
JopNlProp3DAcoVTIDenQ_DEO2_FDTD(; kwargs...)
JopLnProp3DAcoVTIDenQ_DEO2_FDTD(; m₀, kwargs...)

Create a Jets nonlinear or linearized operator for 3D pseudo- visco-acoustic, vertical transverse isotropy, variable density modeling.

Model Parameters

This propagator operates with four model parameters, as shown in the table below.

ParameterDescription
vP wave velocity
ϵModified Thomsen's weak anisotropy parameter
ηModified Alkhalifah's weak anisotropy parameter η
bbuoyancy (reciprocal density)

Pseudo-acoustic approximation

With the pseudo-acoustic approximation we employ, the transformation to self-adjoint form requires a parameter f representing the average ratio of shear-wave to p-wave velocity, and a modfication to Alkhalifah's weak anisotropy parameter η. We show formulas for these two parameters below. The parameter f is specified as a scalar for the entire model, with default value 0.85, implying a shear velocity around 38% of the p wave velocity.

  • η = sqrt[2(ϵ-δ)/(f+2ϵ)]
  • f = 1 - Vₛ² / Vₚ²

For more information about the pseudo-acoustic approximation we employ please see https://library.seg.org/doi/10.1190/segam2016-13878451.1.

Active and passive model parameters

There are two modes of operation that define different sets of active and passive parameters, velocity only and velocity and anisotropy.

An active parameter can be inverted for using the Jacobian linearized operator machinery, and a passive parameter is constant. The two modes are:

ModeActive ParametersPassive Parameters
velocity only[v][ϵ, η, b]
velocity and anisotropy[v, ϵ, η][b]

To make a parameter passive, you pass the array for that parameter in the constructor for the operator, implying that it is in state and does not change with the action of the operator. Parameters that are not passed to the constructor are assumed to be active, and will be part of the model that the operator acts on, stored as a 3D array with the following indices: [Z coord][X coord][Active Parameter]

Examples

Model and acquisition geometry setup

  1. load modules Jets, WaveFD, and JetPackWaveFD
  2. set up the model discretization, coordinate size and spacing
  3. set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
  4. create constant models for v, ϵ, η, b
using Jets, WaveFD, JetPackWaveFD
nz,ny,nx = 100, 80, 60                        # spatial discretization size
dz,dy,dx = 20.0, 20.0, 20.0                   # spatial discretization sampling
ntrec = 1101                                  # number of temporal samples in recorded data
dtrec = 0.0100                                # sample rate for recorded data
dtmod = 0.0025                                # sample rate for modeled data
wavelet = WaveletCausalRicker(f=5.0)          # type of wavelet to use (source signature) 
sz = dz                                       # source Z location 
sy = dy*(ny/2)                                # source Y location 
sx = dx*(nx/2)                                # source X location 
rz = [dz for iy = 1:ny, ix = 1:nx][:];        # Array of receiver Y locations 
ry = [(iy-1)*dy for iy = 1:ny, ix = 1:nx][:]; # Array of receiver Y locations
rx = [(ix-1)*dx for iy = 1:ny, ix=1:nx][:];   # Array of receiver X locations
b = ones(Float32, nz, ny, nx);                # constant buoyancy
ϵ = 0.1 * ones(Float32, nz, ny, nx);          # constant epsilon
η = 0.2 * ones(Float32, nz, ny, nx);          # constant eta

Construct and apply the nonlinear operator (v is active; ϵ, η, b are passive)

  1. create the nonlinear operator F
  2. create the constant velocity model m₀
  3. perform nonlinear forward modeling with constant velocity model m₀ and return the resulting modeled data in d
F = JopNlProp3DAcoVTIDenQ_DEO2_FDTD(; b = b, ϵ = ϵ, η = η, isinterior=true, nsponge = 10, 
    ntrec = ntrec, dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, 
    wavelet = wavelet, sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,:,1] .= 1500;
d = F*m₀;              # forward nonlinear op

Construct and apply the nonlinear operator (v, ϵ, η are active; b is passive)

  1. create the nonlinear operator F
  2. create the model vector and set parameters v, ϵ, η
  3. perform nonlinear forward modeling with model m₀ and return the resulting modeled data in d
F = JopNlProp3DAcoVTIDenQ_DEO2_FDTD(; b = b, isinterior=true, nsponge = 10, 
    ntrec = ntrec, dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, 
    wavelet = wavelet, sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,:,1] .= 1500;
m₀[:,:,:,2] .= 0.1;
m₀[:,:,:,3] .= 0.2;
d = F*m₀;              # forward nonlinear op

Construct and apply the linearized Jacobian operator (method 1)

For this example we assume in the model that v, ϵ, η are active parameters, and b is passive.

  1. create the nonlinear operator F directly by construction of a jet.
  2. create the model vector and set parameters v, ϵ, η
  3. create the Jacobian operator J by linearization of F at point m₀
  4. create a random model perturbation vector δm.
  5. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  6. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm.

Note that the Jacobian operators J and J' require the serialized nonlinear forward wavefield, and this is generated automatically whenever required. If you watch the logging to standard out from this example, you will first see the finite difference evolution for the nonlinear forward, followed by the linearized forward, and finally the linearized adjoint.

F = JopNlProp3DAcoVTIDenQ_DEO2_FDTD(; b = b, isinterior=true, nsponge = 10, ntrec = ntrec, 
    dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, 
    sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,:,1] .= 1500;
m₀[:,:,:,2] .= 0.1;
m₀[:,:,:,3] .= 0.2;
J = jacobian(F, m₀)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Construct and apply the linearized Jacobian operator (method 2)

For this example we assume in the model that v, ϵ, η are active parameters, and b is passive.

  1. create the constant velocity model m₀
  2. create the Jacobian operator J at point m₀ directly by construction of a jet.
  3. create a random model perturbation vector δm.
  4. perform linearized forward (Born) modeling on the model perturbation vector δm and return the resulting data perturbation in δd.
  5. perform linearized adjoint (Born) migration on the data perturbation vector δd and return the resulting model perturbation in δm. 
m₀ = zeros(Float32, nz, ny, nx, 3);
m₀[:,:,:,1] .= 1500;
m₀[:,:,:,2] .= 0.1;
m₀[:,:,:,3] .= 0.2;
J = JopLnProp3DAcoVTIDenQ_DEO2_FDTD(; m₀ = m₀, b = b, isinterior=true, nsponge = 10, 
    ntrec = ntrec, dz = dz, dy = dy, dx = dx, dtrec = dtrec, dtmod = dtmod, 
    wavelet = wavelet, sz = sz, sy = sy, sx = sx, ry = ry, rz = rz, rx = rx)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op

Required Parameters

  • m₀ the point at which the jet is linearized. Note this argument is required in the constuctor for JopLnProp2DAcoVTIDenQ_DEO2_FDTD but not JopNlProp2DAcoVTIDenQ_DEO2_FDTD. This constuctor is shown in the last example above. Please note that you must consider which parameters are active and passive, per the discussion on model parameters above and examples.
  • dtmod The sample rate for the modeled data. You can establish a lower bound for the modeling sample rate with the following expression: dt = 0.75 * 0.38 * max(dx,dz) / maximum(v). Note that the usual Courant–Friedrichs–Lewy condition (CFL condition) for stability in finite difference modeling is modified heuristically to include the impact of the visco-acoustic implementation of this operator, requiring a 25% smaller smaple rate.
  • dtrec The number of time samples in the recorded data. Note requirement that dtrec be an even multiple of dtmod.
  • ntrec The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined from ntrec, dtrec, and dtmod.

Optional Parameters

Defaults for arguments are shown inside square brackets.

  • b [ones(Float32,0,0,0)] the buoyancy (reciprocal density) array.
  • f [0.85] See the discussion above regarding the Pseudo-acoustic approximation used here
  • srcfieldfile [joinpath(tempdir(), "field-e9ebc6be-af0c-4c5a-b7dc-5b5ccb6774a4.bin")] the full path to a scratch file used for the serializationof the compressed nonlinear source wavefield.
  • comptype [nothing] the type of compression to use for the serialization of the nonlinear source wavefield. The type of compression must be one of:
    • nothing - no compression.
    • Float32 - if wavefield is Float64, then do a simple conversion to Float32.
    • UInt32 - compression using CvxCompress (windowing + 2D wavelet transform + thresholding + quantization + run-length-encoding).
  • compscale [1e-2] determines the thresholding for the compression of the nonlinear source wavefield prior to serialization. Larger values mean more aggressive compression. You can likely increase compscale to 1.0 before you start to notice differences in output from Jacobian operations.
  • nz_subcube, ny_subcube, nx_subcube [32] The Z, Y, and X sizes of windows used for compression of the nonlinear source wavefield with CvxCompress. Note the requirement [8 <= n*_subcube <=256].
  • isinterior [false] boolean flag that indicates how the nonlinear source wavefield is saved. For large models, operation will be faster with isinterior = true, but the linearization correctness test may fail.
    • true the entire model including absorbing boundaries is serialized and deserialized
    • false the interior part of the model excluding absorbing boundaries is serialized and deserialized
  • sz [0.0] Array of source Z coordinate. Note that if multiple sources are provided, they will will be injected simultaneously during finite difference evolution.
  • sy [0.0] Array of source Y coordinate.
  • sx [0.0] Array of source X coordinate.
  • st [0.0] Array of source delay times.
  • interpmethod [:hicks] Type of physical interpolation for sources and receivers. For locations that are not on the physical grid coordinates, interpolation is used either to inject or extract information. interpmethod must be one of:
    • :hicks Hicks 3D sinc interpolation (up to 8x8x8 nonzero points per location)
    • :linear bilinear interpolation (up to 2x2x2 nonzero points per location)
  • rz [[0.0]] 2D array of receiver Z coordinates
  • ry [[0.0]] 2D array of receiver Y coordinates
  • rx [[0.0]] 2D array of receiver Z coordinates
  • z0 [0.0] Origin of physical coordinates for the Z dimension.
  • y0 [0.0] Origin of physical coordinates for the Y dimension.
  • x0 [0.0] Origin of physical coordinates for the X dimension.
  • dz [10.0] Spacing of physical coordinates in the Z dimension.
  • dy [10.0] Spacing of physical coordinates in the Y dimension.
  • dx [10.0] Spacing of physical coordinates in the X dimension.
  • freqQ [5.0] The center frequency for the Maxwell body approximation to dissipation only attenuation. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qMin [0.1] The minimum value for Qp at the boundary of the model used in our Maxwell body approximation to dissipation only attenuation. This is not a physically meaningful value for Qp, as we use the attenuation to implement absorbing boundary conditions and eliminate outgoing waves on the boundaries of the computational domain. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • qInterior [100.0] the value for Qp in the interior of the model used in our Maxwell body approximation to dissipation only attenuation. This is the value for Qp away from the absorbing boundaries and is a physically meaningful value. Please see JetPackWaveFD package documentation for more information concerning the attenuation model.
  • padz, pady, padx [0.0], [0.0], [0.0] - apply extra padding to the survey determined aperture in Ginsu. Please see Ginsu for more information.
  • nbz_cache, nby_cache, nbx_cache [512], [8], [8] The size of cache blocks in the Z, X, and Y dimensions. In general the cache block in the Z (fast) dimension should be ≥ the entire size of that dimension, and the cache block size in the slower dimensions is generally small in order to allow the entire block to fit in cache.
  • nsponge [50] The number of grid cells to use for the absorbing boundary. For high fidelity modeling this should be > 60 grid points, but can be significantly smaller for some use cases like low frequency full waveform inversion.
  • wavelet [WaveletCausalRicker(f=5.0)] The source wavelet, can be specified as either a Wavelet type or an array.
  • freesurface [false] Determines if a free surface (true) or absorbing (false) top boundary condition is applied.
  • imgcondition ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for RTM). "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint only exists for "standard" imaging condition currently.
  • RTM_weight [0.5] determines the balance of short wavelengths and long wavelengths in the imaging condition. A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition.
  • nthreads [Sys.CPU_THREADS] The number of threads to use for OpenMP parallelization of the modeling.
  • reportinterval [500] The interval at which information about the propagtion is logged.

See also: Ginsu, WaveletSine, WaveletRicker, WaveletMinPhaseRicker, WaveletDerivRicker, WaveletCausalRicker, WaveletOrmsby, WaveletMinPhaseOrmsby

source
JetPackWaveFD.interior — Method
I = interior(ginsu)

Get the index range corresponding to the interior region which excludes the sponge.

source
JetPackWaveFD.interior — Method
mi = interior(ginsu, mg)

Get a view of the Ginsu'd subset mg corresponding to its interior region which excludes the sponge.

source
JetPackWaveFD.srcillum! — Method
srcillum(J)
srcillum(J, m)
srcillum!(y, J, m)

Compute and return the source illumination for Jets operator J. The source illumination is defined as the sum of squares of the wavefield amplitudes everywhere in the model over the time of the finite difference evolution.

If J::Jop is a Jets block operator or distributed block operator, the source illuminations from all blocks will be accumulated with a simple sum.

srcillum(J) creates and returns an array with the source illumination from J::JopLn, using the current location defined in the Jet.

srcillum(J, m) creates and returns an array with the source illumination from J::Jop for the location m.

srcillum!(y, J, m) zeros the passed array y and then accumulates to the source illumination from J::Jop at the location m into y.

source
JetPackWaveFD.srcillum — Method
srcillum(J)
srcillum(J, m)
srcillum!(y, J, m)

Compute and return the source illumination for Jets operator J. The source illumination is defined as the sum of squares of the wavefield amplitudes everywhere in the model over the time of the finite difference evolution.

If J::Jop is a Jets block operator or distributed block operator, the source illuminations from all blocks will be accumulated with a simple sum.

srcillum(J) creates and returns an array with the source illumination from J::JopLn, using the current location defined in the Jet.

srcillum(J, m) creates and returns an array with the source illumination from J::Jop for the location m.

srcillum!(y, J, m) zeros the passed array y and then accumulates to the source illumination from J::Jop at the location m into y.

source
JetPackWaveFD.sub! — Method
sub!(mg, ginsu, m; extend=true, interior=false)

Store the subset of m corresponding to ginsu in mg.

source
JetPackWaveFD.sub — Method
mg = sub(ginsu, m; extend=true, interior=false)

Store the subset of m corresponding to ginsu in mg. New memory is allocated for mg. interior=false includes the sponge region.

source
JetPackWaveFD.super! — Method
super!(m, ginsu, mg; interior=false, accumulate=false)

Store the superset of mg corresponding to ginsu in m. interior=false includes the sponge region.

source
JetPackWaveFD.Ginsu — Method
g = Ginsu(r0, dr, nr, aperturer, ndamp; T=Float32)

Create a Ginsu object from absolute aperture.

required parameters[1,2]

  • aperturer::NTuple{Range} aperture in each dimension
source
JetPackWaveFD.Ginsu — Method
g = Ginsu(r0, dr, nr, sour, recr, padr, ndamp; dims=(:z,:y,:x), stencilhalfwidth=2, vector_width=8)

Create a Ginsu object that defines the model aperture for a given shot. g is used to subset earth models using the sub and sub! methods, and provides the inverse operation using the super, super! and super_accumulate! methods. The padding for the earth model subset is determined from the source and receiver positions, along with the padding padr and damping region ndamp parameters.

required parameters[1,2]

  • r0::NTuple{N,Real} model origin in each model dimension
  • dr::NTuple{N,Real} model cell widths in each model dimension
  • nr::NTuple{N,Int} model cell counts in each model dimension
  • sour::NTuple{N,Vector{Float}} source locations in each model dimension
  • recr::NTuple{N,Vector{Float}} receiver locations in each model dimension
  • padr::NTuple{N,NTuple{Real,Real}} padding beyond the survey aperture in each model dimension[3,4]
  • ndamp::NTuple{N,NTuple{Int,Int}} damping region for absorbing boundaies in each model dimension[3,4]

named optional parameters

  • dim=(:z,:y,:x) axis ordering with the default being z fast, and x slow
  • stencilhalfwidth=0 if there is a free-surface, set stencilhalfwidth, padr, and ndamp accordingly. This will add stencilhalfwidth cells above the free surface to allow copying the mirrored model to implement the free surface boundary condition.
  • vector_width=8 sets the width required for vectorization of the code, and ensures that the size(ginsu) divided by vector_width has remainder zero for each dimension.

type specification

  • g.lextrng is the logical (1-based) indices for the Ginsu model subset
  • g.lintrng is the logical (1-based) indices for the Ginsu model subset interior
  • rₒ is the physical origin of the Ginsu model subset
  • δr is the grid cell sizes

Notes

  1. N is the number of model dimensions
  2. rₒ[i] is the model origin along the ith model dimension
  3. the model padding is the combination of padr and ndamp
  4. padr[i][1] is the padding for the dimension start, and padr[i][2] for the end. The same is true for ndamp.

Example (free-surface)

o---------------c---x--------------x---c-------o
|               |   |              |   |       |
|               |   |              |   |       |
|               |   |              |   |       |
|               |   |              |   |       |
|               |   x--------------x   |       |
|               |                      |       |
o---------------c----------------------c-------o
  • o denotes the original model
  • c denotes the Ginsu module subset
  • x denotes the Ginsu model subset interior
  • the region falling in-between the Ginsu model subset and the Ginsu model subset interior is the absorbing boundary region.
source

Index