From 4a87288bdb1050dcf06c9d5b99b366dee00ea2ff Mon Sep 17 00:00:00 2001 From: MatthewEai Date: Thu, 15 Jun 2023 12:14:20 +0000 Subject: [PATCH 1/3] added functionality for mixed imaging condition --- src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl | 45 ++++++++++++++++++++++++--- src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl | 42 ++++++++++++++++++++++--- src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl | 42 ++++++++++++++++++++++--- src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl | 42 ++++++++++++++++++++++--- src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl | 42 ++++++++++++++++++++++--- src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl | 43 ++++++++++++++++++++++--- 6 files changed, 231 insertions(+), 25 deletions(-) diff --git a/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl index 894a368..6baebae 100644 --- a/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl +++ b/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl @@ -33,6 +33,7 @@ function JetProp2DAcoIsoDenQ_DEO2_FDTD(; wavelet = WaveletCausalRicker(f=5.0), freesurface = false, imgcondition = "standard", + RTM_weight = 0.5, nthreads = Sys.CPU_THREADS, reportinterval = 500) # active and passive earth model properties. The active set is in the model-space @@ -126,10 +127,11 @@ function JetProp2DAcoIsoDenQ_DEO2_FDTD(; icdict = Dict( lowercase("standard") => WaveFD.ImagingConditionStandard(), lowercase("FWI") => WaveFD.ImagingConditionWaveFieldSeparationFWI(), - lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM()) + lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM(), + lowercase("MIX") => WaveFD.ImagingConditionWaveFieldSeparationMIX()) if lowercase(imgcondition) ∉ keys(icdict) - error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]") + error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM, MIX]") end Jet( @@ -170,6 +172,7 @@ function JetProp2DAcoIsoDenQ_DEO2_FDTD(; nbx_cache = nbx_cache, wavelet = wavelet, freesurface = freesurface, + RTM_weight, imgcondition = get(icdict, lowercase(imgcondition), WaveFD.ImagingConditionStandard()), nthreads = nthreads, reportinterval = reportinterval, @@ -356,7 +359,11 @@ Defaults for arguments are shown inside square brackets. * `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). Note the true adjoint only exists for "standard" imaging condition currently. + 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` 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. @@ -722,10 +729,21 @@ function JopProp2DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end δm_ginsu = Dict{String,Array{Float32,2}}() + + # Temporary model updates if mixed imaging condition used. + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM = Dict{String,Array{Float32,2}}() + δm_All = Dict{String,Array{Float32,2}}() + end + for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, nx_ginsu) + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM[prop] = zeros(Float32, nz_ginsu, nx_ginsu) + δm_All[prop] = zeros(Float32, nz_ginsu, nx_ginsu) + end end - + ginsu_interior_range = interior(kwargs[:ginsu]) # Get receiver interpolation coefficients @@ -797,7 +815,24 @@ function JopProp2DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() + # Standard born accumulation for FWI, RTM, or standard imaging condition + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + else + # Born accumulation for mix imaging condition + imcon = WaveFD.ImagingConditionStandard() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) + + imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) + + weightAll = 1.0f0 - kwargs[:RTM_weight] + weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 + + for prop in keys(kwargs[:active_modelset]) + δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) + end + end end end for prop in keys(kwargs[:active_modelset]) diff --git a/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl index d14d95f..afd0162 100644 --- a/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl @@ -35,6 +35,7 @@ function JetProp2DAcoTTIDenQ_DEO2_FDTD(; wavelet = WaveletCausalRicker(f=5.0), freesurface = false, imgcondition = "standard", + RTM_weight = 0.5, nthreads = Sys.CPU_THREADS, reportinterval = 500) # active and passive earth model properties. The active set is in the model-space @@ -116,10 +117,11 @@ function JetProp2DAcoTTIDenQ_DEO2_FDTD(; icdict = Dict( lowercase("standard") => WaveFD.ImagingConditionStandard(), lowercase("FWI") => WaveFD.ImagingConditionWaveFieldSeparationFWI(), - lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM()) + lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM(), + lowercase("MIX") => WaveFD.ImagingConditionWaveFieldSeparationMIX()) if lowercase(imgcondition) ∉ keys(icdict) - error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]") + error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM, MIX]") end # construct: @@ -163,6 +165,7 @@ function JetProp2DAcoTTIDenQ_DEO2_FDTD(; wavelet = wavelet, freesurface = freesurface, imgcondition = get(icdict, lowercase(imgcondition), WaveFD.ImagingConditionStandard()), + RTM_weight = RTM_weight, nthreads = nthreads, reportinterval = reportinterval, stats = Dict{String,Float64}("MCells/s"=>0.0, "%io"=>0.0, "%inject/extract"=>0.0, "%imaging"=>0.0))) @@ -401,7 +404,11 @@ Defaults for arguments are shown inside square brackets. * `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). Note the true adjoint only exists for "standard" imaging condition currently. + 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` 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. @@ -763,8 +770,18 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA earth_ginsu["f"] .= kwargs[:f] δm_ginsu = Dict{String,Array{Float32,2}}() + # Temporary model updates if mixed imaging condition used. + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM = Dict{String,Array{Float32,2}}() + δm_All = Dict{String,Array{Float32,2}}() + end + for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, nx_ginsu) + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM[prop] = zeros(Float32, nz_ginsu, nx_ginsu) + δm_All[prop] = zeros(Float32, nz_ginsu, nx_ginsu) + end end ginsu_interior_range = interior(kwargs[:ginsu]) @@ -835,7 +852,24 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() + # Standard born accumulation for FWI, RTM, or standard imaging condition + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + else + # Born accumulation for mix imaging condition + imcon = WaveFD.ImagingConditionStandard() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) + + imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) + + weightAll = 1.0f0 - kwargs[:RTM_weight] + weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 + + for prop in keys(kwargs[:active_modelset]) + δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) + end + end end end for prop in keys(kwargs[:active_modelset]) diff --git a/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl index 84e9cd3..28246f7 100644 --- a/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl @@ -34,6 +34,7 @@ function JetProp2DAcoVTIDenQ_DEO2_FDTD(; wavelet = WaveletCausalRicker(f=5.0), freesurface = false, imgcondition = "standard", + RTM_weight = 0.5, nthreads = Sys.CPU_THREADS, reportinterval = 500) # active and passive earth model properties. The active set is in the model-space @@ -113,10 +114,11 @@ function JetProp2DAcoVTIDenQ_DEO2_FDTD(; icdict = Dict( lowercase("standard") => WaveFD.ImagingConditionStandard(), lowercase("FWI") => WaveFD.ImagingConditionWaveFieldSeparationFWI(), - lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM()) + lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM(), + lowercase("MIX") => WaveFD.ImagingConditionWaveFieldSeparationMIX()) if lowercase(imgcondition) ∉ keys(icdict) - error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]") + error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM, MIX]") end Jet( @@ -159,6 +161,7 @@ function JetProp2DAcoVTIDenQ_DEO2_FDTD(; wavelet = wavelet, freesurface = freesurface, imgcondition = get(icdict, lowercase(imgcondition), WaveFD.ImagingConditionStandard()), + RTM_weight = RTM_weight, nthreads = nthreads, reportinterval = reportinterval, stats = Dict{String,Float64}("MCells/s"=>0.0, "%io"=>0.0, "%inject/extract"=>0.0, "%imaging"=>0.0))) @@ -393,7 +396,11 @@ Defaults for arguments are shown inside square brackets. * `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). Note the true adjoint only exists for "standard" imaging condition currently. + 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` 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. @@ -762,8 +769,18 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA earth_ginsu["f"] .= kwargs[:f] δm_ginsu = Dict{String,Array{Float32,2}}() + # Temporary model updates if mixed imaging condition used. + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM = Dict{String,Array{Float32,2}}() + δm_All = Dict{String,Array{Float32,2}}() + end + for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, nx_ginsu) + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM[prop] = zeros(Float32, nz_ginsu, nx_ginsu) + δm_All[prop] = zeros(Float32, nz_ginsu, nx_ginsu) + end end ginsu_interior_range = interior(kwargs[:ginsu]) @@ -836,7 +853,24 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() + # Standard born accumulation for FWI, RTM, or standard imaging condition + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + else + # Born accumulation for mix imaging condition + imcon = WaveFD.ImagingConditionStandard() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) + + imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) + + weightAll = 1.0f0 - kwargs[:RTM_weight] + weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 + + for prop in keys(kwargs[:active_modelset]) + δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) + end + end end end for prop in keys(kwargs[:active_modelset]) diff --git a/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl index 929f6bd..309619a 100644 --- a/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl +++ b/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl @@ -41,6 +41,7 @@ function JetProp3DAcoIsoDenQ_DEO2_FDTD(; wavelet = WaveletCausalRicker(f=5.0), freesurface = false, imgcondition = "standard", + RTM_weight = 0.5, nthreads = Sys.CPU_THREADS, reportinterval = 500) @@ -137,10 +138,11 @@ function JetProp3DAcoIsoDenQ_DEO2_FDTD(; icdict = Dict( lowercase("standard") => WaveFD.ImagingConditionStandard(), lowercase("FWI") => WaveFD.ImagingConditionWaveFieldSeparationFWI(), - lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM()) + lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM(), + lowercase("MIX") => WaveFD.ImagingConditionWaveFieldSeparationMIX()) if lowercase(imgcondition) ∉ keys(icdict) - error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]") + error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM, MIX]") end Jet( @@ -187,6 +189,7 @@ function JetProp3DAcoIsoDenQ_DEO2_FDTD(; wavelet = wavelet, freesurface = freesurface, imgcondition = get(icdict, lowercase(imgcondition), WaveFD.ImagingConditionStandard()), + RTM_weight = RTM_weight, nthreads = nthreads, reportinterval = reportinterval, stats = Dict{String,Float64}("MCells/s"=>0.0, "%io"=>0.0, "%inject/extract"=>0.0, "%imaging"=>0.0))) @@ -378,7 +381,11 @@ Defaults for arguments are shown inside square brackets. * `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). Note the true adjoint only exists for "standard" imaging condition currently. + 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` 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. @@ -752,8 +759,18 @@ function JopProp3DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end δm_ginsu = Dict{String,Array{Float32,3}}() + # Temporary model updates if mixed imaging condition used. + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM = Dict{String,Array{Float32,3}}() + δm_All = Dict{String,Array{Float32,3}}() + end + for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + δm_All[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + end end ginsu_interior_range = interior(kwargs[:ginsu]) @@ -829,7 +846,24 @@ function JopProp3DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() + # Standard born accumulation for FWI, RTM, or standard imaging condition + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + else + # Born accumulation for mix imaging condition + imcon = WaveFD.ImagingConditionStandard() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) + + imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) + + weightAll = 1.0f0 - kwargs[:RTM_weight] + weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 + + for prop in keys(kwargs[:active_modelset]) + δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) + end + end end end for prop in keys(kwargs[:active_modelset]) diff --git a/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl index 2a5900e..0d009fa 100644 --- a/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl @@ -43,6 +43,7 @@ function JetProp3DAcoTTIDenQ_DEO2_FDTD(; wavelet = WaveletCausalRicker(f=5.0), freesurface = false, imgcondition = "standard", + RTM_weight = 0.5, nthreads = Sys.CPU_THREADS, reportinterval = 500) @@ -128,10 +129,11 @@ function JetProp3DAcoTTIDenQ_DEO2_FDTD(; icdict = Dict( lowercase("standard") => WaveFD.ImagingConditionStandard(), lowercase("FWI") => WaveFD.ImagingConditionWaveFieldSeparationFWI(), - lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM()) + lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM(), + lowercase("MIX") => WaveFD.ImagingConditionWaveFieldSeparationMIX()) if lowercase(imgcondition) ∉ keys(icdict) - error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]") + error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM, MIX]") end # construct: @@ -180,6 +182,7 @@ function JetProp3DAcoTTIDenQ_DEO2_FDTD(; wavelet = wavelet, freesurface = freesurface, imgcondition = get(icdict, lowercase(imgcondition), WaveFD.ImagingConditionStandard()), + RTM_weight = RTM_weight, nthreads = nthreads, reportinterval = reportinterval, stats = Dict{String,Float64}("MCells/s"=>0.0, "%io"=>0.0, "%inject/extract"=>0.0, "%imaging"=>0.0))) @@ -432,7 +435,11 @@ Defaults for arguments are shown inside square brackets. * `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). Note the true adjoint only exists for "standard" imaging condition currently. + 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` 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. @@ -808,8 +815,18 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA earth_ginsu["f"] .= kwargs[:f] δm_ginsu = Dict{String,Array{Float32,3}}() + # Temporary model updates if mixed imaging condition used. + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM = Dict{String,Array{Float32,3}}() + δm_All = Dict{String,Array{Float32,3}}() + end + for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + δm_All[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + end end ginsu_interior_range = interior(kwargs[:ginsu]) @@ -884,7 +901,24 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() + # Standard born accumulation for FWI, RTM, or standard imaging condition + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + else + # Born accumulation for mix imaging condition + imcon = WaveFD.ImagingConditionStandard() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) + + imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) + + weightAll = 1.0f0 - kwargs[:RTM_weight] + weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 + + for prop in keys(kwargs[:active_modelset]) + δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) + end + end end end for prop in keys(kwargs[:active_modelset]) diff --git a/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl index 223bb5d..76e7a38 100644 --- a/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl @@ -41,6 +41,7 @@ function JetProp3DAcoVTIDenQ_DEO2_FDTD(; wavelet = WaveletCausalRicker(f=5.0), freesurface = false, imgcondition = "standard", + RTM_weight = 0.5, nthreads = Sys.CPU_THREADS, reportinterval = 500) @@ -122,10 +123,11 @@ function JetProp3DAcoVTIDenQ_DEO2_FDTD(; icdict = Dict( lowercase("standard") => WaveFD.ImagingConditionStandard(), lowercase("FWI") => WaveFD.ImagingConditionWaveFieldSeparationFWI(), - lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM()) + lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM(), + lowercase("MIX") => WaveFD.ImagingConditionWaveFieldSeparationMIX()) if lowercase(imgcondition) ∉ keys(icdict) - error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]") + error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM, MIX]") end # construct: @@ -174,6 +176,7 @@ function JetProp3DAcoVTIDenQ_DEO2_FDTD(; wavelet = wavelet, freesurface = freesurface, imgcondition = get(icdict, lowercase(imgcondition), WaveFD.ImagingConditionStandard()), + RTM_weight = RTM_weight, nthreads = nthreads, reportinterval = reportinterval, stats = Dict{String,Float64}("MCells/s"=>0.0, "%io"=>0.0, "%inject/extract"=>0.0, "%imaging"=>0.0))) @@ -418,7 +421,11 @@ Defaults for arguments are shown inside square brackets. * `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). Note the true adjoint only exists for "standard" imaging condition currently. + 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` 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. @@ -791,8 +798,18 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA earth_ginsu["f"] .= kwargs[:f] δm_ginsu = Dict{String,Array{Float32,3}}() + # Temporary model updates if mixed imaging condition used. + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM = Dict{String,Array{Float32,3}}() + δm_All = Dict{String,Array{Float32,3}}() + end + for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + δm_All[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + end end ginsu_interior_range = interior(kwargs[:ginsu]) @@ -867,7 +884,25 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() + # Standard born accumulation for FWI, RTM, or standard imaging condition + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + else + # Born accumulation for mix imaging condition + imcon = WaveFD.ImagingConditionStandard() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) + + imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) + + weightAll = 1.0f0 - kwargs[:RTM_weight] + weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 + + for prop in keys(kwargs[:active_modelset]) + δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) + end + end + end end for prop in keys(kwargs[:active_modelset]) From 900b8354a2ba6916082a3e55144e2e565fe55f52 Mon Sep 17 00:00:00 2001 From: samkaplan Date: Fri, 4 Aug 2023 16:06:38 +0000 Subject: [PATCH 2/3] mixed imaging conditions --- src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl | 42 ++++++++------------------- src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl | 39 +++++++------------------ src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl | 39 +++++++------------------ src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl | 39 +++++++------------------ src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl | 39 +++++++------------------ src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl | 40 +++++++------------------ 6 files changed, 67 insertions(+), 171 deletions(-) diff --git a/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl index 6baebae..4623e41 100644 --- a/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl +++ b/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl @@ -729,21 +729,14 @@ function JopProp2DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end δm_ginsu = Dict{String,Array{Float32,2}}() - - # Temporary model updates if mixed imaging condition used. - if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() - δm_RTM = Dict{String,Array{Float32,2}}() - δm_All = Dict{String,Array{Float32,2}}() - end - for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, nx_ginsu) - if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() - δm_RTM[prop] = zeros(Float32, nz_ginsu, nx_ginsu) - δm_All[prop] = zeros(Float32, nz_ginsu, nx_ginsu) + if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) + δm_ginsu["rtm_$prop"] = zeros(Float32, nz_ginsu, nx_ginsu) + δm_ginsu["all_$prop"] = zeros(Float32, nz_ginsu, nx_ginsu) end end - + ginsu_interior_range = interior(kwargs[:ginsu]) # Get receiver interpolation coefficients @@ -815,30 +808,19 @@ function JopProp2DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() - # Standard born accumulation for FWI, RTM, or standard imaging condition - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) - else - # Born accumulation for mix imaging condition - imcon = WaveFD.ImagingConditionStandard() - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) - - imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) - - weightAll = 1.0f0 - kwargs[:RTM_weight] - weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 - - for prop in keys(kwargs[:active_modelset]) - δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) - end - end + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) end end + set_zero_subnormals(false) + for prop in keys(kwargs[:active_modelset]) + if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) + weight_all = 1 - kwargs[:RTM_weight] + weight_short = 2*kwargs[:RTM_weight] - 1 + δm_ginsu[prop] .= (δm_ginsu["all_$prop"] .* weight_all) .+ (δm_ginsu["rtm_$prop"] .* weight_short) + end δm_ginsu[prop] .*= kwargs[:dtrec] / kwargs[:dtmod] end - set_zero_subnormals(false) JopProp2DAcoIsoDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im) diff --git a/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl index afd0162..b15c323 100644 --- a/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl @@ -770,17 +770,11 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA earth_ginsu["f"] .= kwargs[:f] δm_ginsu = Dict{String,Array{Float32,2}}() - # Temporary model updates if mixed imaging condition used. - if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() - δm_RTM = Dict{String,Array{Float32,2}}() - δm_All = Dict{String,Array{Float32,2}}() - end - for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, nx_ginsu) - if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() - δm_RTM[prop] = zeros(Float32, nz_ginsu, nx_ginsu) - δm_All[prop] = zeros(Float32, nz_ginsu, nx_ginsu) + if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) + δm_ginsu["rtm_$prop"] = zeros(Float32, nz_ginsu, nx_ginsu) + δm_ginsu["all_$prop"] = zeros(Float32, nz_ginsu, nx_ginsu) end end @@ -852,30 +846,19 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() - # Standard born accumulation for FWI, RTM, or standard imaging condition - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) - else - # Born accumulation for mix imaging condition - imcon = WaveFD.ImagingConditionStandard() - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) - - imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) - - weightAll = 1.0f0 - kwargs[:RTM_weight] - weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 - - for prop in keys(kwargs[:active_modelset]) - δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) - end - end + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) end end + set_zero_subnormals(false) + for prop in keys(kwargs[:active_modelset]) + if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) + weight_all = 1 - kwargs[:RTM_weight] + weight_short = 2*kwargs[:RTM_weight] - 1 + δm_ginsu[prop] .= (δm_ginsu["all_$prop"] .* weight_all) .+ (δm_ginsu["rtm_$prop"] .* weight_short) + end δm_ginsu[prop] .*= kwargs[:dtrec] / kwargs[:dtmod] end - set_zero_subnormals(false) JopProp2DAcoTTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im) # undo ginsu diff --git a/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl index 28246f7..83d2119 100644 --- a/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl @@ -769,17 +769,11 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA earth_ginsu["f"] .= kwargs[:f] δm_ginsu = Dict{String,Array{Float32,2}}() - # Temporary model updates if mixed imaging condition used. - if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() - δm_RTM = Dict{String,Array{Float32,2}}() - δm_All = Dict{String,Array{Float32,2}}() - end - for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, nx_ginsu) - if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() - δm_RTM[prop] = zeros(Float32, nz_ginsu, nx_ginsu) - δm_All[prop] = zeros(Float32, nz_ginsu, nx_ginsu) + if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) + δm_ginsu["rtm_$prop"] = zeros(Float32, nz_ginsu, nx_ginsu) + δm_ginsu["all_$prop"] = zeros(Float32, nz_ginsu, nx_ginsu) end end @@ -853,30 +847,19 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() - # Standard born accumulation for FWI, RTM, or standard imaging condition - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) - else - # Born accumulation for mix imaging condition - imcon = WaveFD.ImagingConditionStandard() - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) - - imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) - - weightAll = 1.0f0 - kwargs[:RTM_weight] - weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 - - for prop in keys(kwargs[:active_modelset]) - δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) - end - end + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) end end + set_zero_subnormals(false) + for prop in keys(kwargs[:active_modelset]) + if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) + weight_all = 1 - kwargs[:RTM_weight] + weight_short = 2*kwargs[:RTM_weight] - 1 + δm_ginsu[prop] .= (δm_ginsu["all_$prop"] .* weight_all) .+ (δm_ginsu["rtm_$prop"] .* weight_short) + end δm_ginsu[prop] .*= kwargs[:dtrec] / kwargs[:dtmod] end - set_zero_subnormals(false) JopProp2DAcoVTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im) # undo ginsu diff --git a/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl index 309619a..bf5a7aa 100644 --- a/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl +++ b/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl @@ -759,17 +759,11 @@ function JopProp3DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end δm_ginsu = Dict{String,Array{Float32,3}}() - # Temporary model updates if mixed imaging condition used. - if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() - δm_RTM = Dict{String,Array{Float32,3}}() - δm_All = Dict{String,Array{Float32,3}}() - end - for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) - if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() - δm_RTM[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) - δm_All[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) + δm_ginsu["rtm_$prop"] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + δm_ginsu["all_$prop"] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) end end @@ -846,30 +840,19 @@ function JopProp3DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() - # Standard born accumulation for FWI, RTM, or standard imaging condition - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) - else - # Born accumulation for mix imaging condition - imcon = WaveFD.ImagingConditionStandard() - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) - - imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) - - weightAll = 1.0f0 - kwargs[:RTM_weight] - weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 - - for prop in keys(kwargs[:active_modelset]) - δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) - end - end + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) end end + set_zero_subnormals(false) + for prop in keys(kwargs[:active_modelset]) + if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) + weight_all = 1 - kwargs[:RTM_weight] + weight_short = 2*kwargs[:RTM_weight] - 1 + δm_ginsu[prop] .= (δm_ginsu["all_$prop"] .* weight_all) .+ (δm_ginsu["rtm_$prop"] .* weight_short) + end δm_ginsu[prop] .*= kwargs[:dtrec] / kwargs[:dtmod] end - set_zero_subnormals(false) JopProp3DAcoIsoDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im) diff --git a/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl index 0d009fa..605886a 100644 --- a/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl @@ -815,17 +815,11 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA earth_ginsu["f"] .= kwargs[:f] δm_ginsu = Dict{String,Array{Float32,3}}() - # Temporary model updates if mixed imaging condition used. - if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() - δm_RTM = Dict{String,Array{Float32,3}}() - δm_All = Dict{String,Array{Float32,3}}() - end - for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) - if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() - δm_RTM[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) - δm_All[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) + δm_ginsu["rtm_$prop"] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + δm_ginsu["all_$prop"] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) end end @@ -901,30 +895,19 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() - # Standard born accumulation for FWI, RTM, or standard imaging condition - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) - else - # Born accumulation for mix imaging condition - imcon = WaveFD.ImagingConditionStandard() - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) - - imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) - - weightAll = 1.0f0 - kwargs[:RTM_weight] - weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 - - for prop in keys(kwargs[:active_modelset]) - δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) - end - end + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) end end + set_zero_subnormals(false) + for prop in keys(kwargs[:active_modelset]) + if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) + weight_all = 1 - kwargs[:RTM_weight] + weight_short = 2*kwargs[:RTM_weight] - 1 + δm_ginsu[prop] .= (δm_ginsu["all_$prop"] .* weight_all) .+ (δm_ginsu["rtm_$prop"] .* weight_short) + end δm_ginsu[prop] .*= kwargs[:dtrec] / kwargs[:dtmod] end - set_zero_subnormals(false) JopProp3DAcoTTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im) # undo ginsu diff --git a/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl index 76e7a38..fcd2537 100644 --- a/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl @@ -798,17 +798,11 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA earth_ginsu["f"] .= kwargs[:f] δm_ginsu = Dict{String,Array{Float32,3}}() - # Temporary model updates if mixed imaging condition used. - if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() - δm_RTM = Dict{String,Array{Float32,3}}() - δm_All = Dict{String,Array{Float32,3}}() - end - for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) - if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() - δm_RTM[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) - δm_All[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) + δm_ginsu["rtm_$prop"] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + δm_ginsu["all_$prop"] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) end end @@ -884,31 +878,19 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() - # Standard born accumulation for FWI, RTM, or standard imaging condition - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) - else - # Born accumulation for mix imaging condition - imcon = WaveFD.ImagingConditionStandard() - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) - - imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) - - weightAll = 1.0f0 - kwargs[:RTM_weight] - weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 - - for prop in keys(kwargs[:active_modelset]) - δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) - end - end - + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) end end + set_zero_subnormals(false) + for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] .*= kwargs[:dtrec] / kwargs[:dtmod] + if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) + weight_all = 1 - kwargs[:RTM_weight] + weight_short = 2*kwargs[:RTM_weight] - 1 + δm_ginsu[prop] .= (δm_ginsu["all_$prop"] .* weight_all) .+ (δm_ginsu["rtm_$prop"] .* weight_short) + end end - set_zero_subnormals(false) JopProp3DAcoVTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im) # undo ginsu From a4742ee8f328b9d79baab74a7477a1bebd306715 Mon Sep 17 00:00:00 2001 From: samkaplan Date: Fri, 4 Aug 2023 16:06:54 +0000 Subject: [PATCH 3/3] Revert "mixed imaging conditions" This reverts commit 900b8354a2ba6916082a3e55144e2e565fe55f52. --- src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl | 42 +++++++++++++++++++-------- src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl | 39 ++++++++++++++++++------- src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl | 39 ++++++++++++++++++------- src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl | 39 ++++++++++++++++++------- src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl | 39 ++++++++++++++++++------- src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl | 40 ++++++++++++++++++------- 6 files changed, 171 insertions(+), 67 deletions(-) diff --git a/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl index 4623e41..6baebae 100644 --- a/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl +++ b/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl @@ -729,14 +729,21 @@ function JopProp2DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end δm_ginsu = Dict{String,Array{Float32,2}}() + + # Temporary model updates if mixed imaging condition used. + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM = Dict{String,Array{Float32,2}}() + δm_All = Dict{String,Array{Float32,2}}() + end + for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, nx_ginsu) - if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) - δm_ginsu["rtm_$prop"] = zeros(Float32, nz_ginsu, nx_ginsu) - δm_ginsu["all_$prop"] = zeros(Float32, nz_ginsu, nx_ginsu) + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM[prop] = zeros(Float32, nz_ginsu, nx_ginsu) + δm_All[prop] = zeros(Float32, nz_ginsu, nx_ginsu) end end - + ginsu_interior_range = interior(kwargs[:ginsu]) # Get receiver interpolation coefficients @@ -808,19 +815,30 @@ function JopProp2DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() + # Standard born accumulation for FWI, RTM, or standard imaging condition + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + else + # Born accumulation for mix imaging condition + imcon = WaveFD.ImagingConditionStandard() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) + + imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) + + weightAll = 1.0f0 - kwargs[:RTM_weight] + weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 + + for prop in keys(kwargs[:active_modelset]) + δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) + end + end end end - set_zero_subnormals(false) - for prop in keys(kwargs[:active_modelset]) - if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) - weight_all = 1 - kwargs[:RTM_weight] - weight_short = 2*kwargs[:RTM_weight] - 1 - δm_ginsu[prop] .= (δm_ginsu["all_$prop"] .* weight_all) .+ (δm_ginsu["rtm_$prop"] .* weight_short) - end δm_ginsu[prop] .*= kwargs[:dtrec] / kwargs[:dtmod] end + set_zero_subnormals(false) JopProp2DAcoIsoDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im) diff --git a/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl index b15c323..afd0162 100644 --- a/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl @@ -770,11 +770,17 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA earth_ginsu["f"] .= kwargs[:f] δm_ginsu = Dict{String,Array{Float32,2}}() + # Temporary model updates if mixed imaging condition used. + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM = Dict{String,Array{Float32,2}}() + δm_All = Dict{String,Array{Float32,2}}() + end + for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, nx_ginsu) - if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) - δm_ginsu["rtm_$prop"] = zeros(Float32, nz_ginsu, nx_ginsu) - δm_ginsu["all_$prop"] = zeros(Float32, nz_ginsu, nx_ginsu) + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM[prop] = zeros(Float32, nz_ginsu, nx_ginsu) + δm_All[prop] = zeros(Float32, nz_ginsu, nx_ginsu) end end @@ -846,19 +852,30 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() + # Standard born accumulation for FWI, RTM, or standard imaging condition + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + else + # Born accumulation for mix imaging condition + imcon = WaveFD.ImagingConditionStandard() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) + + imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) + + weightAll = 1.0f0 - kwargs[:RTM_weight] + weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 + + for prop in keys(kwargs[:active_modelset]) + δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) + end + end end end - set_zero_subnormals(false) - for prop in keys(kwargs[:active_modelset]) - if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) - weight_all = 1 - kwargs[:RTM_weight] - weight_short = 2*kwargs[:RTM_weight] - 1 - δm_ginsu[prop] .= (δm_ginsu["all_$prop"] .* weight_all) .+ (δm_ginsu["rtm_$prop"] .* weight_short) - end δm_ginsu[prop] .*= kwargs[:dtrec] / kwargs[:dtmod] end + set_zero_subnormals(false) JopProp2DAcoTTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im) # undo ginsu diff --git a/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl index 83d2119..28246f7 100644 --- a/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl @@ -769,11 +769,17 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA earth_ginsu["f"] .= kwargs[:f] δm_ginsu = Dict{String,Array{Float32,2}}() + # Temporary model updates if mixed imaging condition used. + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM = Dict{String,Array{Float32,2}}() + δm_All = Dict{String,Array{Float32,2}}() + end + for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, nx_ginsu) - if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) - δm_ginsu["rtm_$prop"] = zeros(Float32, nz_ginsu, nx_ginsu) - δm_ginsu["all_$prop"] = zeros(Float32, nz_ginsu, nx_ginsu) + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM[prop] = zeros(Float32, nz_ginsu, nx_ginsu) + δm_All[prop] = zeros(Float32, nz_ginsu, nx_ginsu) end end @@ -847,19 +853,30 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() + # Standard born accumulation for FWI, RTM, or standard imaging condition + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + else + # Born accumulation for mix imaging condition + imcon = WaveFD.ImagingConditionStandard() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) + + imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) + + weightAll = 1.0f0 - kwargs[:RTM_weight] + weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 + + for prop in keys(kwargs[:active_modelset]) + δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) + end + end end end - set_zero_subnormals(false) - for prop in keys(kwargs[:active_modelset]) - if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) - weight_all = 1 - kwargs[:RTM_weight] - weight_short = 2*kwargs[:RTM_weight] - 1 - δm_ginsu[prop] .= (δm_ginsu["all_$prop"] .* weight_all) .+ (δm_ginsu["rtm_$prop"] .* weight_short) - end δm_ginsu[prop] .*= kwargs[:dtrec] / kwargs[:dtmod] end + set_zero_subnormals(false) JopProp2DAcoVTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im) # undo ginsu diff --git a/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl index bf5a7aa..309619a 100644 --- a/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl +++ b/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl @@ -759,11 +759,17 @@ function JopProp3DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end δm_ginsu = Dict{String,Array{Float32,3}}() + # Temporary model updates if mixed imaging condition used. + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM = Dict{String,Array{Float32,3}}() + δm_All = Dict{String,Array{Float32,3}}() + end + for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) - if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) - δm_ginsu["rtm_$prop"] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) - δm_ginsu["all_$prop"] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + δm_All[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) end end @@ -840,19 +846,30 @@ function JopProp3DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() + # Standard born accumulation for FWI, RTM, or standard imaging condition + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + else + # Born accumulation for mix imaging condition + imcon = WaveFD.ImagingConditionStandard() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) + + imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) + + weightAll = 1.0f0 - kwargs[:RTM_weight] + weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 + + for prop in keys(kwargs[:active_modelset]) + δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) + end + end end end - set_zero_subnormals(false) - for prop in keys(kwargs[:active_modelset]) - if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) - weight_all = 1 - kwargs[:RTM_weight] - weight_short = 2*kwargs[:RTM_weight] - 1 - δm_ginsu[prop] .= (δm_ginsu["all_$prop"] .* weight_all) .+ (δm_ginsu["rtm_$prop"] .* weight_short) - end δm_ginsu[prop] .*= kwargs[:dtrec] / kwargs[:dtmod] end + set_zero_subnormals(false) JopProp3DAcoIsoDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im) diff --git a/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl index 605886a..0d009fa 100644 --- a/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl @@ -815,11 +815,17 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA earth_ginsu["f"] .= kwargs[:f] δm_ginsu = Dict{String,Array{Float32,3}}() + # Temporary model updates if mixed imaging condition used. + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM = Dict{String,Array{Float32,3}}() + δm_All = Dict{String,Array{Float32,3}}() + end + for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) - if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) - δm_ginsu["rtm_$prop"] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) - δm_ginsu["all_$prop"] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + δm_All[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) end end @@ -895,19 +901,30 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() + # Standard born accumulation for FWI, RTM, or standard imaging condition + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + else + # Born accumulation for mix imaging condition + imcon = WaveFD.ImagingConditionStandard() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) + + imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) + + weightAll = 1.0f0 - kwargs[:RTM_weight] + weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 + + for prop in keys(kwargs[:active_modelset]) + δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) + end + end end end - set_zero_subnormals(false) - for prop in keys(kwargs[:active_modelset]) - if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) - weight_all = 1 - kwargs[:RTM_weight] - weight_short = 2*kwargs[:RTM_weight] - 1 - δm_ginsu[prop] .= (δm_ginsu["all_$prop"] .* weight_all) .+ (δm_ginsu["rtm_$prop"] .* weight_short) - end δm_ginsu[prop] .*= kwargs[:dtrec] / kwargs[:dtmod] end + set_zero_subnormals(false) JopProp3DAcoTTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im) # undo ginsu diff --git a/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl index fcd2537..76e7a38 100644 --- a/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl @@ -798,11 +798,17 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA earth_ginsu["f"] .= kwargs[:f] δm_ginsu = Dict{String,Array{Float32,3}}() + # Temporary model updates if mixed imaging condition used. + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM = Dict{String,Array{Float32,3}}() + δm_All = Dict{String,Array{Float32,3}}() + end + for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) - if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) - δm_ginsu["rtm_$prop"] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) - δm_ginsu["all_$prop"] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + if kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationMIX() + δm_RTM[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) + δm_All[prop] = zeros(Float32, nz_ginsu, ny_ginsu, nx_ginsu) end end @@ -878,19 +884,31 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA end # born accumulation - cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX() + # Standard born accumulation for FWI, RTM, or standard imaging condition + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields) + else + # Born accumulation for mix imaging condition + imcon = WaveFD.ImagingConditionStandard() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields) + + imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM() + cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields) + + weightAll = 1.0f0 - kwargs[:RTM_weight] + weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0 + + for prop in keys(kwargs[:active_modelset]) + δm_ginsu[prop] = (δm_All[prop] .* weightAll) .+ (δm_RTM[prop] .* weightShort) + end + end + end end - set_zero_subnormals(false) - for prop in keys(kwargs[:active_modelset]) δm_ginsu[prop] .*= kwargs[:dtrec] / kwargs[:dtmod] - if isa(kwargs[:imgcondition], WaveFD.ImagingConditionWaveFieldSeparationMIX) - weight_all = 1 - kwargs[:RTM_weight] - weight_short = 2*kwargs[:RTM_weight] - 1 - δm_ginsu[prop] .= (δm_ginsu["all_$prop"] .* weight_all) .+ (δm_ginsu["rtm_$prop"] .* weight_short) - end end + set_zero_subnormals(false) JopProp3DAcoVTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im) # undo ginsu