Reference
JetPackWaveFD.JopLnProp2DAcoIsoDenQ_DEO2_FDTD
— MethodJopNlProp2DAcoIsoDenQ_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.
Parameter | Description |
---|---|
v | P wave velocity |
b | buoyancy (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
- load modules Jets, WaveFD, and JetPackWaveFD
- set up the model discretization, coordinate size and spacing
- set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
- 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
- create the nonlinear operator
F
- create the constant velocity model m₀
- perform nonlinear forward modeling with constant velocity model
m₀
and return the resulting modeled data ind
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)
- create the nonlinear operator
F
directly by construction of a jet. - create the constant velocity model m₀
- create the Jacobian operator
J
by linearization ofF
at pointm₀
- create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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)
- create the constant velocity model m₀
- create the Jacobian operator
J
at pointm₀
directly by construction of a jet. - create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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 forJopLnProp2DAcoIsoDenQ_DEO2_FDTD
but notJopNlProp2DAcoIsoDenQ_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 thatdtrec
be an even multiple ofdtmod
.ntrec
The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined fromntrec
,dtrec
, anddtmod
.
Optional Parameters
Defaults for arguments are shown inside square brackets.
b [ones(Float32,0,0)]
the buoyancy (reciprocal density) array.srcfieldfile [joinpath(tempdir(), "field-1d507251-ccaa-4ead-8ab3-f918b4d8dde6.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 isFloat64
, 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 withCvxCompress
. 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 withisinterior = true
, but the linearization correctness test may fail.true
the entire model including absorbing boundaries is serialized and deserializedfalse
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 coordinatesrx [[0.0]]
2D array of receiver Z coordinatesz0 [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 seeJetPackWaveFD
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 seeJetPackWaveFD
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 seeJetPackWaveFD
package documentation for more information concerning the attenuation model.padz, padx [0.0], [0.0]
- apply extra padding to the survey determined aperture inGinsu
. Please seeGinsu
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
JetPackWaveFD.JopLnProp2DAcoTTIDenQ_DEO2_FDTD
— MethodJopNlProp2DAcoTTIDenQ_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.
Parameter | Description |
---|---|
v | P wave velocity |
ϵ | Modified Thomsen's weak anisotropy parameter |
η | Modified Alkhalifah's weak anisotropy parameter η |
b | buoyancy (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:
Mode | Active Parameters | Passive 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
- load modules Jets, WaveFD, and JetPackWaveFD
- set up the model discretization, coordinate size and spacing
- set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
- 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)
- create the nonlinear operator
F
- create the model vector and set parameter
v
- perform nonlinear forward modeling with model
m₀
and return the resulting modeled data ind
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)
- create the nonlinear operator
F
- create the model vector and set parameters
v, ϵ, η
- perform nonlinear forward modeling with model
m₀
and return the resulting modeled data ind
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.
- create the nonlinear operator
F
directly by construction of a jet. - create the model vector and set parameters
v, ϵ, η
- create the Jacobian operator
J
by linearization ofF
at pointm₀
- create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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.
- create the model vector and set parameters
v, ϵ, η
- create the Jacobian operator
J
at pointm₀
directly by construction of a jet. - create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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 forJopLnProp2DAcoTTIDenQ_DEO2_FDTD
but notJopNlProp2DAcoTTIDenQ_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 thatdtrec
be an even multiple ofdtmod
.ntrec
The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined fromntrec
,dtrec
, anddtmod
.
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 buoyancyb
.f [0.85]
See the discussion above regarding the Pseudo-acoustic approximation used heresrcfieldfile [joinpath(tempdir(), "field-e01027f8-fafa-4ae5-a510-c04225ac1607.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 isFloat64
, 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 withCvxCompress
. 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 withisinterior = true
, but the linearization correctness test may fail.true
the entire model including absorbing boundaries is serialized and deserializedfalse
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 coordinatesrx [[0.0]]
2D array of receiver Z coordinatesz0 [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 seeJetPackWaveFD
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 seeJetPackWaveFD
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 seeJetPackWaveFD
package documentation for more information concerning the attenuation model.padz, padx [0.0], [0.0]
- apply extra padding to the survey determined aperture inGinsu
. Please seeGinsu
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
JetPackWaveFD.JopLnProp2DAcoVTIDenQ_DEO2_FDTD
— MethodJopNlProp2DAcoVTIDenQ_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.
Parameter | Description |
---|---|
v | P wave velocity |
ϵ | Modified Thomsen's weak anisotropy parameter |
η | Modified Alkhalifah's weak anisotropy parameter η |
b | buoyancy (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:
Mode | Active Parameters | Passive 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
- load modules Jets, WaveFD, and JetPackWaveFD
- set up the model discretization, coordinate size and spacing
- set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
- 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)
- create the nonlinear operator
F
- create the model vector and set parameter
v
- perform nonlinear forward modeling with model
m₀
and return the resulting modeled data ind
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)
- create the nonlinear operator
F
- create the model vector and set parameters
v, ϵ, η
- perform nonlinear forward modeling with model
m₀
and return the resulting modeled data ind
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.
- create the nonlinear operator
F
directly by construction of a jet. - create the model vector and set parameters
v, ϵ, η
- create the Jacobian operator
J
by linearization ofF
at pointm₀
- create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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.
- create the model vector and set parameters
v, ϵ, η
- create the Jacobian operator
J
at pointm₀
directly by construction of a jet. - create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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 forJopLnProp2DAcoVTIDenQ_DEO2_FDTD
but notJopNlProp2DAcoVTIDenQ_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 thatdtrec
be an even multiple ofdtmod
.ntrec
The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined fromntrec
,dtrec
, anddtmod
.
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 heresrcfieldfile [joinpath(tempdir(), "field-7ff0fa67-17d0-4de9-a7da-ff8210be0156.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 isFloat64
, 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 withCvxCompress
. 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 withisinterior = true
, but the linearization correctness test may fail.true
the entire model including absorbing boundaries is serialized and deserializedfalse
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 coordinatesrx [[0.0]]
2D array of receiver Z coordinatesz0 [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 seeJetPackWaveFD
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 seeJetPackWaveFD
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 seeJetPackWaveFD
package documentation for more information concerning the attenuation model.padz, padx [0.0], [0.0]
- apply extra padding to the survey determined aperture inGinsu
. Please seeGinsu
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
JetPackWaveFD.JopLnProp3DAcoIsoDenQ_DEO2_FDTD
— MethodJopNlProp3DAcoIsoDenQ_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.
Parameter | Description |
---|---|
v | P wave velocity |
b | buoyancy (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
- load modules Jets, WaveFD, and JetPackWaveFD
- set up the model discretization, coordinate size and spacing
- set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
- 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
- create the nonlinear operator
F
- create the constant velocity model m₀
- perform nonlinear forward modeling with constant velocity model
m₀
and return the resulting modeled data ind
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)
- create the nonlinear operator
F
directly by construction of a jet. - create the constant velocity model m₀
- create the Jacobian operator
J
by linearization ofF
at pointm₀
- create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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)
- create the constant velocity model m₀
- create the Jacobian operator
J
at pointm₀
directly by construction of a jet. - create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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 forJopLnProp3DAcoIsoDenQ_DEO2_FDTD
but notJopNlProp3DAcoIsoDenQ_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 thatdtrec
be an even multiple ofdtmod
.ntrec
The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined fromntrec
,dtrec
, anddtmod
.
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-4f642538-9ee8-4172-afe4-9d58f5b6de53.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 isFloat64
, 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 withCvxCompress
. 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 withisinterior = true
, but the linearization correctness test may fail.true
the entire model including absorbing boundaries is serialized and deserializedfalse
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 coordinatesry [[0.0]]
2D array of receiver Y coordinatesrx [[0.0]]
2D array of receiver Z coordinatesz0 [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 seeJetPackWaveFD
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 seeJetPackWaveFD
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 seeJetPackWaveFD
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 inGinsu
. Please seeGinsu
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
JetPackWaveFD.JopLnProp3DAcoTTIDenQ_DEO2_FDTD
— MethodJopNlProp3DAcoTTIDenQ_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.
Parameter | Description |
---|---|
v | P wave velocity |
ϵ | Modified Thomsen's weak anisotropy parameter |
η | Modified Alkhalifah's weak anisotropy parameter η |
b | buoyancy (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:
Mode | Active Parameters | Passive 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
- load modules Jets, WaveFD, and JetPackWaveFDFDFD
- set up the model discretization, coordinate size and spacing
- set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
- 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)
- create the nonlinear operator
F
- create the constant velocity model m₀
- perform nonlinear forward modeling with constant velocity model
m₀
and return the resulting modeled data ind
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)
- create the nonlinear operator
F
- create the model vector and set parameters
v, ϵ, η
- perform nonlinear forward modeling with model
m₀
and return the resulting modeled data ind
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.
- create the nonlinear operator
F
directly by construction of a jet. - create the model vector and set parameters
v, ϵ, η
- create the Jacobian operator
J
by linearization ofF
at pointm₀
- create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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.
- create the constant velocity model m₀
- create the Jacobian operator
J
at pointm₀
directly by construction of a jet. - create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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 forJopLnProp2DAcoTTIDenQ_DEO2_FDTD
but notJopNlProp2DAcoTTIDenQ_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 thatdtrec
be an even multiple ofdtmod
.ntrec
The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined fromntrec
,dtrec
, anddtmod
.
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 buoyancyb
.f [0.85]
See the discussion above regarding the Pseudo-acoustic approximation used heresrcfieldfile ["field-0849ed16-3e1b-4c10-be0b-c177cf6366ff.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 isFloat64
, 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 withCvxCompress
. 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 withisinterior = true
, but the linearization correctness test may fail.true
the entire model including absorbing boundaries is serialized and deserializedfalse
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 coordinatesry [[0.0]]
2D array of receiver Y coordinatesrx [[0.0]]
2D array of receiver Z coordinatesz0 [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 seeJetPackWaveFDFDFD
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 seeJetPackWaveFDFDFD
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 seeJetPackWaveFDFDFD
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 inGinsu
. Please seeGinsu
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
JetPackWaveFD.JopLnProp3DAcoVTIDenQ_DEO2_FDTD
— MethodJopNlProp3DAcoVTIDenQ_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.
Parameter | Description |
---|---|
v | P wave velocity |
ϵ | Modified Thomsen's weak anisotropy parameter |
η | Modified Alkhalifah's weak anisotropy parameter η |
b | buoyancy (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:
Mode | Active Parameters | Passive 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
- load modules Jets, WaveFD, and JetPackWaveFD
- set up the model discretization, coordinate size and spacing
- set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
- 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)
- create the nonlinear operator
F
- create the constant velocity model m₀
- perform nonlinear forward modeling with constant velocity model
m₀
and return the resulting modeled data ind
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)
- create the nonlinear operator
F
- create the model vector and set parameters
v, ϵ, η
- perform nonlinear forward modeling with model
m₀
and return the resulting modeled data ind
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.
- create the nonlinear operator
F
directly by construction of a jet. - create the model vector and set parameters
v, ϵ, η
- create the Jacobian operator
J
by linearization ofF
at pointm₀
- create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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.
- create the constant velocity model m₀
- create the Jacobian operator
J
at pointm₀
directly by construction of a jet. - create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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 forJopLnProp2DAcoVTIDenQ_DEO2_FDTD
but notJopNlProp2DAcoVTIDenQ_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 thatdtrec
be an even multiple ofdtmod
.ntrec
The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined fromntrec
,dtrec
, anddtmod
.
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 heresrcfieldfile [joinpath(tempdir(), "field-b130513e-9cfe-45bc-98a9-6657bf0785b1.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 isFloat64
, 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 withCvxCompress
. 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 withisinterior = true
, but the linearization correctness test may fail.true
the entire model including absorbing boundaries is serialized and deserializedfalse
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 coordinatesry [[0.0]]
2D array of receiver Y coordinatesrx [[0.0]]
2D array of receiver Z coordinatesz0 [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 seeJetPackWaveFD
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 seeJetPackWaveFD
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 seeJetPackWaveFD
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 inGinsu
. Please seeGinsu
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
JetPackWaveFD.JopNlProp2DAcoIsoDenQ_DEO2_FDTD
— MethodJopNlProp2DAcoIsoDenQ_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.
Parameter | Description |
---|---|
v | P wave velocity |
b | buoyancy (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
- load modules Jets, WaveFD, and JetPackWaveFD
- set up the model discretization, coordinate size and spacing
- set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
- 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
- create the nonlinear operator
F
- create the constant velocity model m₀
- perform nonlinear forward modeling with constant velocity model
m₀
and return the resulting modeled data ind
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)
- create the nonlinear operator
F
directly by construction of a jet. - create the constant velocity model m₀
- create the Jacobian operator
J
by linearization ofF
at pointm₀
- create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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)
- create the constant velocity model m₀
- create the Jacobian operator
J
at pointm₀
directly by construction of a jet. - create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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 forJopLnProp2DAcoIsoDenQ_DEO2_FDTD
but notJopNlProp2DAcoIsoDenQ_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 thatdtrec
be an even multiple ofdtmod
.ntrec
The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined fromntrec
,dtrec
, anddtmod
.
Optional Parameters
Defaults for arguments are shown inside square brackets.
b [ones(Float32,0,0)]
the buoyancy (reciprocal density) array.srcfieldfile [joinpath(tempdir(), "field-1d507251-ccaa-4ead-8ab3-f918b4d8dde6.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 isFloat64
, 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 withCvxCompress
. 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 withisinterior = true
, but the linearization correctness test may fail.true
the entire model including absorbing boundaries is serialized and deserializedfalse
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 coordinatesrx [[0.0]]
2D array of receiver Z coordinatesz0 [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 seeJetPackWaveFD
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 seeJetPackWaveFD
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 seeJetPackWaveFD
package documentation for more information concerning the attenuation model.padz, padx [0.0], [0.0]
- apply extra padding to the survey determined aperture inGinsu
. Please seeGinsu
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
JetPackWaveFD.JopNlProp2DAcoTTIDenQ_DEO2_FDTD
— MethodJopNlProp2DAcoTTIDenQ_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.
Parameter | Description |
---|---|
v | P wave velocity |
ϵ | Modified Thomsen's weak anisotropy parameter |
η | Modified Alkhalifah's weak anisotropy parameter η |
b | buoyancy (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:
Mode | Active Parameters | Passive 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
- load modules Jets, WaveFD, and JetPackWaveFD
- set up the model discretization, coordinate size and spacing
- set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
- 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)
- create the nonlinear operator
F
- create the model vector and set parameter
v
- perform nonlinear forward modeling with model
m₀
and return the resulting modeled data ind
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)
- create the nonlinear operator
F
- create the model vector and set parameters
v, ϵ, η
- perform nonlinear forward modeling with model
m₀
and return the resulting modeled data ind
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.
- create the nonlinear operator
F
directly by construction of a jet. - create the model vector and set parameters
v, ϵ, η
- create the Jacobian operator
J
by linearization ofF
at pointm₀
- create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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.
- create the model vector and set parameters
v, ϵ, η
- create the Jacobian operator
J
at pointm₀
directly by construction of a jet. - create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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 forJopLnProp2DAcoTTIDenQ_DEO2_FDTD
but notJopNlProp2DAcoTTIDenQ_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 thatdtrec
be an even multiple ofdtmod
.ntrec
The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined fromntrec
,dtrec
, anddtmod
.
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 buoyancyb
.f [0.85]
See the discussion above regarding the Pseudo-acoustic approximation used heresrcfieldfile [joinpath(tempdir(), "field-e01027f8-fafa-4ae5-a510-c04225ac1607.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 isFloat64
, 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 withCvxCompress
. 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 withisinterior = true
, but the linearization correctness test may fail.true
the entire model including absorbing boundaries is serialized and deserializedfalse
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 coordinatesrx [[0.0]]
2D array of receiver Z coordinatesz0 [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 seeJetPackWaveFD
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 seeJetPackWaveFD
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 seeJetPackWaveFD
package documentation for more information concerning the attenuation model.padz, padx [0.0], [0.0]
- apply extra padding to the survey determined aperture inGinsu
. Please seeGinsu
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
JetPackWaveFD.JopNlProp2DAcoVTIDenQ_DEO2_FDTD
— MethodJopNlProp2DAcoVTIDenQ_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.
Parameter | Description |
---|---|
v | P wave velocity |
ϵ | Modified Thomsen's weak anisotropy parameter |
η | Modified Alkhalifah's weak anisotropy parameter η |
b | buoyancy (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:
Mode | Active Parameters | Passive 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
- load modules Jets, WaveFD, and JetPackWaveFD
- set up the model discretization, coordinate size and spacing
- set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
- 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)
- create the nonlinear operator
F
- create the model vector and set parameter
v
- perform nonlinear forward modeling with model
m₀
and return the resulting modeled data ind
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)
- create the nonlinear operator
F
- create the model vector and set parameters
v, ϵ, η
- perform nonlinear forward modeling with model
m₀
and return the resulting modeled data ind
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.
- create the nonlinear operator
F
directly by construction of a jet. - create the model vector and set parameters
v, ϵ, η
- create the Jacobian operator
J
by linearization ofF
at pointm₀
- create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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.
- create the model vector and set parameters
v, ϵ, η
- create the Jacobian operator
J
at pointm₀
directly by construction of a jet. - create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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 forJopLnProp2DAcoVTIDenQ_DEO2_FDTD
but notJopNlProp2DAcoVTIDenQ_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 thatdtrec
be an even multiple ofdtmod
.ntrec
The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined fromntrec
,dtrec
, anddtmod
.
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 heresrcfieldfile [joinpath(tempdir(), "field-7ff0fa67-17d0-4de9-a7da-ff8210be0156.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 isFloat64
, 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 withCvxCompress
. 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 withisinterior = true
, but the linearization correctness test may fail.true
the entire model including absorbing boundaries is serialized and deserializedfalse
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 coordinatesrx [[0.0]]
2D array of receiver Z coordinatesz0 [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 seeJetPackWaveFD
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 seeJetPackWaveFD
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 seeJetPackWaveFD
package documentation for more information concerning the attenuation model.padz, padx [0.0], [0.0]
- apply extra padding to the survey determined aperture inGinsu
. Please seeGinsu
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
JetPackWaveFD.JopNlProp3DAcoIsoDenQ_DEO2_FDTD
— MethodJopNlProp3DAcoIsoDenQ_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.
Parameter | Description |
---|---|
v | P wave velocity |
b | buoyancy (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
- load modules Jets, WaveFD, and JetPackWaveFD
- set up the model discretization, coordinate size and spacing
- set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
- 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
- create the nonlinear operator
F
- create the constant velocity model m₀
- perform nonlinear forward modeling with constant velocity model
m₀
and return the resulting modeled data ind
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)
- create the nonlinear operator
F
directly by construction of a jet. - create the constant velocity model m₀
- create the Jacobian operator
J
by linearization ofF
at pointm₀
- create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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)
- create the constant velocity model m₀
- create the Jacobian operator
J
at pointm₀
directly by construction of a jet. - create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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 forJopLnProp3DAcoIsoDenQ_DEO2_FDTD
but notJopNlProp3DAcoIsoDenQ_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 thatdtrec
be an even multiple ofdtmod
.ntrec
The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined fromntrec
,dtrec
, anddtmod
.
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-4f642538-9ee8-4172-afe4-9d58f5b6de53.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 isFloat64
, 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 withCvxCompress
. 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 withisinterior = true
, but the linearization correctness test may fail.true
the entire model including absorbing boundaries is serialized and deserializedfalse
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 coordinatesry [[0.0]]
2D array of receiver Y coordinatesrx [[0.0]]
2D array of receiver Z coordinatesz0 [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 seeJetPackWaveFD
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 seeJetPackWaveFD
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 seeJetPackWaveFD
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 inGinsu
. Please seeGinsu
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
JetPackWaveFD.JopNlProp3DAcoTTIDenQ_DEO2_FDTD
— MethodJopNlProp3DAcoTTIDenQ_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.
Parameter | Description |
---|---|
v | P wave velocity |
ϵ | Modified Thomsen's weak anisotropy parameter |
η | Modified Alkhalifah's weak anisotropy parameter η |
b | buoyancy (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:
Mode | Active Parameters | Passive 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
- load modules Jets, WaveFD, and JetPackWaveFDFDFD
- set up the model discretization, coordinate size and spacing
- set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
- 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)
- create the nonlinear operator
F
- create the constant velocity model m₀
- perform nonlinear forward modeling with constant velocity model
m₀
and return the resulting modeled data ind
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)
- create the nonlinear operator
F
- create the model vector and set parameters
v, ϵ, η
- perform nonlinear forward modeling with model
m₀
and return the resulting modeled data ind
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.
- create the nonlinear operator
F
directly by construction of a jet. - create the model vector and set parameters
v, ϵ, η
- create the Jacobian operator
J
by linearization ofF
at pointm₀
- create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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.
- create the constant velocity model m₀
- create the Jacobian operator
J
at pointm₀
directly by construction of a jet. - create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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 forJopLnProp2DAcoTTIDenQ_DEO2_FDTD
but notJopNlProp2DAcoTTIDenQ_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 thatdtrec
be an even multiple ofdtmod
.ntrec
The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined fromntrec
,dtrec
, anddtmod
.
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 buoyancyb
.f [0.85]
See the discussion above regarding the Pseudo-acoustic approximation used heresrcfieldfile ["field-0849ed16-3e1b-4c10-be0b-c177cf6366ff.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 isFloat64
, 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 withCvxCompress
. 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 withisinterior = true
, but the linearization correctness test may fail.true
the entire model including absorbing boundaries is serialized and deserializedfalse
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 coordinatesry [[0.0]]
2D array of receiver Y coordinatesrx [[0.0]]
2D array of receiver Z coordinatesz0 [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 seeJetPackWaveFDFDFD
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 seeJetPackWaveFDFDFD
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 seeJetPackWaveFDFDFD
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 inGinsu
. Please seeGinsu
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
JetPackWaveFD.JopNlProp3DAcoVTIDenQ_DEO2_FDTD
— MethodJopNlProp3DAcoVTIDenQ_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.
Parameter | Description |
---|---|
v | P wave velocity |
ϵ | Modified Thomsen's weak anisotropy parameter |
η | Modified Alkhalifah's weak anisotropy parameter η |
b | buoyancy (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:
Mode | Active Parameters | Passive 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
- load modules Jets, WaveFD, and JetPackWaveFD
- set up the model discretization, coordinate size and spacing
- set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
- 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)
- create the nonlinear operator
F
- create the constant velocity model m₀
- perform nonlinear forward modeling with constant velocity model
m₀
and return the resulting modeled data ind
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)
- create the nonlinear operator
F
- create the model vector and set parameters
v, ϵ, η
- perform nonlinear forward modeling with model
m₀
and return the resulting modeled data ind
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.
- create the nonlinear operator
F
directly by construction of a jet. - create the model vector and set parameters
v, ϵ, η
- create the Jacobian operator
J
by linearization ofF
at pointm₀
- create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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.
- create the constant velocity model m₀
- create the Jacobian operator
J
at pointm₀
directly by construction of a jet. - create a random model perturbation vector
δm
. - perform linearized forward (Born) modeling on the model perturbation vector
δm
and return the resulting data perturbation inδd
. - 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 forJopLnProp2DAcoVTIDenQ_DEO2_FDTD
but notJopNlProp2DAcoVTIDenQ_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 thatdtrec
be an even multiple ofdtmod
.ntrec
The number of time samples in the recorded data. Note that the number of samples in the modeled data is determined fromntrec
,dtrec
, anddtmod
.
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 heresrcfieldfile [joinpath(tempdir(), "field-b130513e-9cfe-45bc-98a9-6657bf0785b1.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 isFloat64
, 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 withCvxCompress
. 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 withisinterior = true
, but the linearization correctness test may fail.true
the entire model including absorbing boundaries is serialized and deserializedfalse
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 coordinatesry [[0.0]]
2D array of receiver Y coordinatesrx [[0.0]]
2D array of receiver Z coordinatesz0 [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 seeJetPackWaveFD
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 seeJetPackWaveFD
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 seeJetPackWaveFD
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 inGinsu
. Please seeGinsu
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
JetPackWaveFD.interior
— MethodI = interior(ginsu)
Get the index range corresponding to the interior region which excludes the sponge.
JetPackWaveFD.interior
— Methodmi = interior(ginsu, mg)
Get a view of the Ginsu'd subset mg
corresponding to its interior region which excludes the sponge.
JetPackWaveFD.srcillum!
— Methodsrcillum(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
.
JetPackWaveFD.srcillum
— Methodsrcillum(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
.
JetPackWaveFD.sub!
— Methodsub!(mg, ginsu, m; extend=true, interior=false)
Store the subset of m
corresponding to ginsu
in mg
.
JetPackWaveFD.sub
— Methodmg = 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.
JetPackWaveFD.super!
— Methodsuper!(m, ginsu, mg; interior=false, accumulate=false)
Store the superset of mg
corresponding to ginsu
in m
. interior=false
includes the sponge region.
JetPackWaveFD.Ginsu
— Methodg = 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
JetPackWaveFD.Ginsu
— Methodg = 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 dimensiondr::NTuple{N,Real}
model cell widths in each model dimensionnr::NTuple{N,Int}
model cell counts in each model dimensionsour::NTuple{N,Vector{Float}}
source locations in each model dimensionrecr::NTuple{N,Vector{Float}}
receiver locations in each model dimensionpadr::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 beingz
fast, andx
slowstencilhalfwidth=0
if there is a free-surface, setstencilhalfwidth
, padr, and ndamp accordingly. This will addstencilhalfwidth
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 subsetg.lintrng
is the logical (1-based) indices for the Ginsu model subset interiorrₒ
is the physical origin of the Ginsu model subsetδr
is the grid cell sizes
Notes
N
is the number of model dimensionsrₒ[i]
is the model origin along the ith model dimension- the model padding is the combination of
padr
andndamp
padr[i][1]
is the padding for the dimension start, andpadr[i][2]
for the end. The same is true forndamp
.
Example (free-surface)
o---------------c---x--------------x---c-------o
| | | | | |
| | | | | |
| | | | | |
| | | | | |
| | x--------------x | |
| | | |
o---------------c----------------------c-------o
o
denotes the original modelc
denotes the Ginsu module subsetx
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.
Index
JetPackWaveFD.Ginsu
JetPackWaveFD.Ginsu
JetPackWaveFD.JopLnProp2DAcoIsoDenQ_DEO2_FDTD
JetPackWaveFD.JopLnProp2DAcoTTIDenQ_DEO2_FDTD
JetPackWaveFD.JopLnProp2DAcoVTIDenQ_DEO2_FDTD
JetPackWaveFD.JopLnProp3DAcoIsoDenQ_DEO2_FDTD
JetPackWaveFD.JopLnProp3DAcoTTIDenQ_DEO2_FDTD
JetPackWaveFD.JopLnProp3DAcoVTIDenQ_DEO2_FDTD
JetPackWaveFD.JopNlProp2DAcoIsoDenQ_DEO2_FDTD
JetPackWaveFD.JopNlProp2DAcoTTIDenQ_DEO2_FDTD
JetPackWaveFD.JopNlProp2DAcoVTIDenQ_DEO2_FDTD
JetPackWaveFD.JopNlProp3DAcoIsoDenQ_DEO2_FDTD
JetPackWaveFD.JopNlProp3DAcoTTIDenQ_DEO2_FDTD
JetPackWaveFD.JopNlProp3DAcoVTIDenQ_DEO2_FDTD
JetPackWaveFD.interior
JetPackWaveFD.interior
JetPackWaveFD.srcillum
JetPackWaveFD.srcillum!
JetPackWaveFD.sub
JetPackWaveFD.sub!
JetPackWaveFD.super!