From e5e1e81abcee0adb522c5bb9a731183f889a6e59 Mon Sep 17 00:00:00 2001 From: samkaplan Date: Wed, 16 Aug 2023 13:40:10 -0500 Subject: [PATCH] add a time taper for the src illum integration --- benchmark/benchmarks.jl | 2 + src/illumination.jl | 75 ++++++++++++++++++-------- src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl | 4 +- src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl | 4 +- src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl | 4 +- src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl | 4 +- src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl | 4 +- src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl | 4 +- test/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl | 5 ++ test/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl | 5 ++ test/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl | 5 ++ test/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl | 5 ++ test/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl | 5 ++ test/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl | 5 ++ 14 files changed, 96 insertions(+), 35 deletions(-) diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index 36d5964..6d2488f 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -262,4 +262,6 @@ for nthreads in _nthreads SUITE["3DAcoTTIDenQ_DEO2_FDTD, adjoint, compress"]["$nthreads threads"] = @benchmarkable mul!(δm,J',δd) setup=(F=f3dtti($nthreads,$(n_3D.z),$(n_3D.y),$(n_3D.x),$(nb_3D.z),$(nb_3D.y),$(nb_3D.x),true,true); m=1500*ones(domain(F)); d=F*m; J=jacobian!(F,m); δm=zeros(domain(J)); δm[div(n_3D.z,2),div(n_3D.y,2),div(n_3D.x,2)] = 100; δd=J*δm) teardown=(close(F)) end +SUITE["srcillum_helper"] = @benchmarkable JetPackWaveFD.srcillum_helper(x, 1.0f0) setup=(x = rand(Float32,700,600,600)) + SUITE \ No newline at end of file diff --git a/src/illumination.jl b/src/illumination.jl index 6f1c384..6fb1928 100644 --- a/src/illumination.jl +++ b/src/illumination.jl @@ -1,44 +1,70 @@ # Jets/WaveFD bridge +struct IllumOnesVector{T} <: AbstractArray{T,1} + n::Int + IllumOnesVector(::Type{T}, n) where {T} = new{T}(n) +end +IllumOnesVector(::Type{T}) where {T} = IllumOnesVector(T, typemax(Int)) + +Base.size(x::IllumOnesVector) = (x.n,) +Base.IndexStyle(::Type{<:IllumOnesVector}) = IndexLinear() +Base.getindex(_::IllumOnesVector{T}, i::Int) where {T} = one(T) +Base.maximum(_::IllumOnesVector{T}) where {T} = one(T) +Base.minimum(_::IllumOnesVector{T}) where {T} = one(T) +Base.extrema(_::IllumOnesVector{T}) where {T} = (one(T),one(T)) + # 2D source illumination function srcillum(filename::AbstractString, compressor::WaveFD.Compressor, interior::Bool, ::Type{T}, ginsu::Ginsu, nz::Int, nx::Int, ntrec::Int, nthreads::Int) where {T} - srcillum!(zeros(T, nz, nx), filename, compressor, ginsu, ntrec, nthreads) + srcillum!(zeros(T, nz, nx), filename, compressor, IllumOnesVector(T, ntrec), interior, ginsu, ntrec, nthreads, mask) end -function srcillum!(γ::AbstractArray{T,2}, filename::AbstractString, compressor::WaveFD.Compressor, interior::Bool, ginsu::Ginsu, ntrec::Int, nthreads::Int) where {T} +function srcillum(filename::AbstractString, compressor::WaveFD.Compressor, time_mask::AbstractVector{T}, interior::Bool, ::Type{T}, ginsu::Ginsu, nz::Int, nx::Int, ntrec::Int, nthreads::Int) where {T} + srcillum!(zeros(T, nz, nx), filename, compressor, time_mask, interior, ginsu, ntrec, nthreads, mask) +end + +function srcillum!(γ::AbstractArray{T,2}, filename::AbstractString, compressor::WaveFD.Compressor, time_mask::AbstractVector{T}, interior::Bool, ginsu::Ginsu, ntrec::Int, nthreads::Int) where {T} io = open(filename) field_ginsu = zeros(T, size(ginsu, interior=interior)) for it = 1:ntrec WaveFD.compressedread!(io, compressor, it, field_ginsu) - WaveFD.srcillum_helper(field_ginsu, nthreads) + srcillum_helper(field_ginsu, time_mask[it]) super!(γ, ginsu, field_ginsu, accumulate=true, interior=interior) end close(io) γ end +function srcillum_helper(field_ginsu, time_mask) + Threads.@threads :static for i in eachindex(field_ginsu) + field_ginsu[i] *= field_ginsu[i] * time_mask + end +end + # 3D source illumination function srcillum(filename::AbstractString, compressor::WaveFD.Compressor, interior::Bool, ::Type{T}, ginsu::Ginsu, nz::Int, ny::Int, nx::Int, ntrec::Int, nthreads::Int) where {T} - srcillum!(zeros(T, nz, ny, nx), filename, compressor, interior, ginsu, ntrec, nthreads) + srcillum!(zeros(T, nz, ny, nx), filename, compressor, IllumOnesVector(T, ntrec), interior, ginsu, ntrec, nthreads) +end + +function srcillum(filename::AbstractString, compressor::WaveFD.Compressor, time_mask::AbstractVector{T}, interior::Bool, ::Type{T}, ginsu::Ginsu, nz::Int, ny::Int, nx::Int, ntrec::Int, nthreads::Int) where {T} + srcillum!(zeros(T, nz, ny, nx), filename, compressor, time_mask, interior, ginsu, ntrec, nthreads) end -function srcillum!(γ::AbstractArray{T,3}, filename::AbstractString, compressor::WaveFD.Compressor, interior::Bool, ginsu::Ginsu, ntrec::Int, nthreads::Int) where {T} +function srcillum!(γ::AbstractArray{T,3}, filename::AbstractString, compressor::WaveFD.Compressor, time_mask::AbstractVector{T}, interior::Bool, ginsu::Ginsu, ntrec::Int, nthreads::Int) where {T} io = open(filename) field_ginsu = zeros(T, size(ginsu, interior=interior)) for it = 1:ntrec WaveFD.compressedread!(io, compressor, it, field_ginsu) - WaveFD.srcillum_helper(field_ginsu, nthreads) + srcillum_helper(field_ginsu, time_mask[it]) super!(γ, ginsu, field_ginsu, accumulate=true, interior=interior) end close(io) γ end - @doc """ - srcillum(J) - srcillum(J, m) - srcillum!(y, J, m) + srcillum(J[; time_mask]) + srcillum(J, m[; time_mask]) + srcillum!(y, J, m[; time_mask]) 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 @@ -55,25 +81,27 @@ 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`. + +A `time_mask` can be used to weight the integration over time. """ -function srcillum(J::JopLn) +function srcillum(J::JopLn{<:Jet{<:JetAbstractSpace{T}}}; time_mask=IllumOnesVector(T)) where {T} s = zeros(eltype(J), size(domain(J))[1:end-1]) - srcillum!(s, J) + srcillum!(s, J; time_mask) end -function srcillum(J::Jop, m::AbstractArray) +function srcillum(J::Jop{<:Jet{<:JetAbstractSpace{T}}}, m::AbstractArray{T}; time_mask=IllumOnesVector(T)) where {T} s = zeros(eltype(J), size(domain(J))[1:end-1]) - srcillum!(s, J, m) + srcillum!(s, J, m; time_mask) end @doc (@doc srcillum(J)) -srcillum!(γ::AbstractArray, J::Jop) = srcillum!(γ, J, jet(J).mₒ) +srcillum!(γ::AbstractArray, J::Jop{<:Jet{<:JetAbstractSpace{T}}}; time_mask=IllumOnesVector(T)) where {T} = srcillum!(γ, J, jet(J).mₒ; time_mask) # composite operators, TODO: the try-catch here is a bit of a kludge -function srcillum!(γ::AbstractArray, A::T, m::AbstractArray) where {D,R,J<:Jet{D,R,typeof(Jets.JetComposite_f!)},T<:Jop{J}} +function srcillum!(γ::AbstractArray, A::Jet{<:JetAbstractSpace{T},<:JetAbstractSpace,typeof(Jets.JetComposite_f!)}, m::AbstractArray; time_mask=IllumOnesVector(T)) where {T} for Aᵢ in state(A).ops try - srcillum!(γ, Aᵢ, m) + srcillum!(γ, Aᵢ, m; time_mask) break catch e if length(m) == 0 @@ -85,27 +113,28 @@ function srcillum!(γ::AbstractArray, A::T, m::AbstractArray) where {D,R,J<:Jet{ end # block operators -function srcillum!(γ::AbstractArray, A::T, m::AbstractArray) where {D,R,J<:Jet{D,R,typeof(Jets.JetBlock_f!)}, T<:Jop{J}} +function srcillum!(γ::AbstractArray, A::Jet{<:JetAbstractSpace{T},<:JetAbstractSpace,typeof(Jets.JetBlock_f!)}, m::AbstractArray; time_mask=IllumOnesVector(T)) where {T} γ .= 0 for Aᵢ in state(A).ops if length(m) == 0 - srcillum!(γ, Aᵢ, jet(Aᵢ).mₒ) + srcillum!(γ, Aᵢ, jet(Aᵢ).mₒ; time_mask) else - srcillum!(γ, Aᵢ, m) + srcillum!(γ, Aᵢ, m; time_mask) end end γ end # distributed block operators -function srcillum!(γ::AbstractArray, A::T, m::AbstractArray) where {D,R,J<:Jet{D,R,typeof(DistributedJets.JetDBlock_f!)},T<:Jop{J}} +function srcillum!(γ::AbstractArray, A::Jet{<:JetAbstractSpace{T},<:JetAbstractSpace,typeof(DistributedJets.JetDBlock_f!)}, m::AbstractArray; time_mask=IllumOnesVector(T)) where {T} γ .= 0 pids = procs(A) γ_partials = ArrayFutures(γ, DistributedJets.addmasterpid(pids)) + time_masks = bcast(time_mask, DistributedJets.addmasterpid(pids)) - _srcillum!(γ_partials, A, m) = begin srcillum!(localpart(γ_partials), localpart(A), m); nothing end + _srcillum!(γ_partials, A, m, time_masks) = begin srcillum!(localpart(γ_partials), localpart(A), m; time_mask=localpart(time_masks)); nothing end @sync for pid in pids - @async remotecall_fetch(_srcillum!, pid, γ_partials, A, m) + @async remotecall_fetch(_srcillum!, pid, γ_partials, A, m, time_tapers) end reduce!(γ_partials) end diff --git a/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl index a0b2b56..0cafc7e 100644 --- a/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl +++ b/src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl @@ -841,7 +841,7 @@ end modelindex(F::Jop{T}, key) where {D,R,T<:Jet{D,R,typeof(JopProp2DAcoIsoDenQ_DEO2_FDTD_f!)}} = state(F).active_modelset[key] -function srcillum!(γ, A::Jop{T}, m::AbstractArray{Float32}) where {D,R,T<:Jet{D,R,typeof(JopProp2DAcoIsoDenQ_DEO2_FDTD_f!)}} +function srcillum!(γ, A::Jop{T}, m::AbstractArray{Float32}; time_mask = IllumOnesVector(Float32)) where {D,R,T<:Jet{D,R,typeof(JopProp2DAcoIsoDenQ_DEO2_FDTD_f!)}} s = state(A) isvalid, _chksum = isvalid_srcfieldfile(jet(A).mₒ, s.srcfieldhost[], s.srcfieldfile*"-pold", s.chksum[]) if !isvalid @@ -850,7 +850,7 @@ function srcillum!(γ, A::Jop{T}, m::AbstractArray{Float32}) where {D,R,T<:Jet{D s.srcfieldhost[] = gethostname() end open(s.compressor["pold"]) - I = srcillum!(γ, s.srcfieldfile*"-pold", s.compressor["pold"], s.isinterior, s.ginsu, s.ntrec, s.nthreads) + I = srcillum!(γ, s.srcfieldfile*"-pold", s.compressor["pold"], time_mask, s.isinterior, s.ginsu, s.ntrec, s.nthreads) close(s.compressor["pold"]) I end diff --git a/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl index 7e4a7c1..77da34d 100644 --- a/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl @@ -878,7 +878,7 @@ end modelindex(F::Jop{T}, key::AbstractString) where {D,R,T<:Jet{D,R,typeof(JopProp2DAcoTTIDenQ_DEO2_FDTD_f!)}} = state(F).active_modelset[key] -function srcillum!(γ, A::T, m::AbstractArray{Float32}) where {D,R,J<:Jet{D,R,typeof(JopProp2DAcoTTIDenQ_DEO2_FDTD_f!)},T<:Jop{J}} +function srcillum!(γ, A::T, m::AbstractArray{Float32}; time_mask = IllumOnesVector(Float32)) where {D,R,J<:Jet{D,R,typeof(JopProp2DAcoTTIDenQ_DEO2_FDTD_f!)},T<:Jop{J}} s = state(A) isvalid, _chksum = isvalid_srcfieldfile(m, s.srcfieldhost[], s.srcfieldfile*"-pold", s.chksum[]) if !isvalid @@ -887,7 +887,7 @@ function srcillum!(γ, A::T, m::AbstractArray{Float32}) where {D,R,J<:Jet{D,R,ty s.srcfieldhost[] = gethostname() end open(s.compressor["pold"]) - I = srcillum!(γ, s.srcfieldfile*"-pold", s.compressor["pold"], s.isinterior, s.ginsu, s.ntrec, s.nthreads) + I = srcillum!(γ, s.srcfieldfile*"-pold", s.compressor["pold"], time_mask, s.isinterior, s.ginsu, s.ntrec, s.nthreads) close(s.compressor["pold"]) I end diff --git a/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl index 307a6d5..09b7727 100644 --- a/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl @@ -879,7 +879,7 @@ end modelindex(F::Jop{T}, key) where {D,R,T<:Jet{D,R,typeof(JopProp2DAcoVTIDenQ_DEO2_FDTD_f!)}} = state(F).active_modelset[key] -function srcillum!(γ, A::Jop{T}, m::AbstractArray{Float32}) where {D,R,T<:Jet{D,R,typeof(JopProp2DAcoVTIDenQ_DEO2_FDTD_f!)}} +function srcillum!(γ, A::Jop{T}, m::AbstractArray{Float32}; time_mask = IllumOnesVector(Float32)) where {D,R,T<:Jet{D,R,typeof(JopProp2DAcoVTIDenQ_DEO2_FDTD_f!)}} s = state(A) isvalid, _chksum = isvalid_srcfieldfile(jet(A).mₒ, s.srcfieldhost[], s.srcfieldfile*"-pold", s.chksum[]) if !isvalid @@ -888,7 +888,7 @@ function srcillum!(γ, A::Jop{T}, m::AbstractArray{Float32}) where {D,R,T<:Jet{D s.srcfieldhost[] = gethostname() end open(s.compressor["pold"]) - I = srcillum!(γ, s.srcfieldfile*"-pold", s.compressor["pold"], s.isinterior, s.ginsu, s.ntrec, s.nthreads) + I = srcillum!(γ, s.srcfieldfile*"-pold", s.compressor["pold"], time_mask, s.isinterior, s.ginsu, s.ntrec, s.nthreads) close(s.compressor["pold"]) I end diff --git a/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl index d8f59cd..544333e 100644 --- a/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl +++ b/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl @@ -873,7 +873,7 @@ end modelindex(F::Jop{T}, key) where {D,R,T<:Jet{D,R,typeof(JopProp3DAcoIsoDenQ_DEO2_FDTD_f!)}} = state(F).active_modelset[key] -function srcillum!(γ, A::Jop{T}, m::AbstractArray{Float32}) where {D,R,T<:Jet{D,R,typeof(JopProp3DAcoIsoDenQ_DEO2_FDTD_f!)}} +function srcillum!(γ, A::Jop{T}, m::AbstractArray{Float32}; time_mask = IllumOnesVector(Float32)) where {D,R,T<:Jet{D,R,typeof(JopProp3DAcoIsoDenQ_DEO2_FDTD_f!)}} s = state(A) isvalid, _chksum = isvalid_srcfieldfile(jet(A).mₒ, s.srcfieldhost[], s.srcfieldfile*"-pold", s.chksum[]) if !isvalid @@ -882,7 +882,7 @@ function srcillum!(γ, A::Jop{T}, m::AbstractArray{Float32}) where {D,R,T<:Jet{D s.srcfieldhost[] = gethostname() end open(s.compressor["pold"]) - I = srcillum!(γ, s.srcfieldfile*"-pold", s.compressor["pold"], s.isinterior, s.ginsu, s.ntrec, s.nthreads) + I = srcillum!(γ, s.srcfieldfile*"-pold", s.compressor["pold"], time_mask, s.isinterior, s.ginsu, s.ntrec, s.nthreads) close(s.compressor["pold"]) I end diff --git a/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl index c52178e..3eb2b80 100644 --- a/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl @@ -927,7 +927,7 @@ end modelindex(F::Jop{T}, key::AbstractString) where {D,R,T<:Jet{D,R,typeof(JopProp3DAcoTTIDenQ_DEO2_FDTD_f!)}} = state(F).active_modelset[key] -function srcillum!(γ, A::T, m::AbstractArray{Float32}) where {D,R,J<:Jet{D,R,typeof(JopProp3DAcoTTIDenQ_DEO2_FDTD_f!)},T<:Jop{J}} +function srcillum!(γ, A::T, m::AbstractArray{Float32}; time_mask = IllumOnesVector(Float32)) where {D,R,J<:Jet{D,R,typeof(JopProp3DAcoTTIDenQ_DEO2_FDTD_f!)},T<:Jop{J}} s = state(A) isvalid, _chksum = isvalid_srcfieldfile(m, s.srcfieldhost[], s.srcfieldfile*"-pold", s.chksum[]) if !isvalid @@ -936,7 +936,7 @@ function srcillum!(γ, A::T, m::AbstractArray{Float32}) where {D,R,J<:Jet{D,R,ty s.srcfieldhost[] = gethostname() end open(s.compressor["pold"]) - I = srcillum!(γ, s.srcfieldfile*"-pold", s.compressor["pold"], s.isinterior, s.ginsu, s.ntrec, s.nthreads) + I = srcillum!(γ, s.srcfieldfile*"-pold", s.compressor["pold"], time_mask, s.isinterior, s.ginsu, s.ntrec, s.nthreads) close(s.compressor["pold"]) I end diff --git a/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl index 57c08ef..98f3223 100644 --- a/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl +++ b/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl @@ -910,7 +910,7 @@ end modelindex(F::Jop{T}, key::AbstractString) where {D,R,T<:Jet{D,R,typeof(JopProp3DAcoVTIDenQ_DEO2_FDTD_f!)}} = state(F).active_modelset[key] -function srcillum!(γ, A::T, m::AbstractArray{Float32}) where {D,R,J<:Jet{D,R,typeof(JopProp3DAcoVTIDenQ_DEO2_FDTD_f!)},T<:Jop{J}} +function srcillum!(γ, A::T, m::AbstractArray{Float32}; time_mask = IllumOnesVector(Float32)) where {D,R,J<:Jet{D,R,typeof(JopProp3DAcoVTIDenQ_DEO2_FDTD_f!)},T<:Jop{J}} s = state(A) isvalid, _chksum = isvalid_srcfieldfile(m, s.srcfieldhost[], s.srcfieldfile*"-pold", s.chksum[]) if !isvalid @@ -919,7 +919,7 @@ function srcillum!(γ, A::T, m::AbstractArray{Float32}) where {D,R,J<:Jet{D,R,ty s.srcfieldhost[] = gethostname() end open(s.compressor["pold"]) - I = srcillum!(γ, s.srcfieldfile*"-pold", s.compressor["pold"], s.isinterior, s.ginsu, s.ntrec, s.nthreads) + I = srcillum!(γ, s.srcfieldfile*"-pold", s.compressor["pold"], time_mask, s.isinterior, s.ginsu, s.ntrec, s.nthreads) close(s.compressor["pold"]) I end diff --git a/test/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl b/test/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl index b083b91..8bfdda8 100644 --- a/test/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl +++ b/test/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl @@ -188,8 +188,13 @@ end J = jacobian(F, m₀); s1 = srcillum(F, m₀); s2 = srcillum(J); + + time_mask = ones(Float32, state(F, :ntrec)) + time_mask[1:div(state(F,:ntrec),2)] .= 0 + s3 = srcillum(J; time_mask) close(F) @test s1 ≈ s2 + @test norm(s3) < norm(s2) @test maximum(s1) > 0 end diff --git a/test/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl b/test/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl index e2096d3..ec4fa38 100644 --- a/test/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl +++ b/test/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl @@ -182,8 +182,13 @@ end J = jacobian(F, m₀); s1 = srcillum(F, m₀); s2 = srcillum(J); + + time_mask = ones(Float32, state(F, :ntrec)) + time_mask[1:div(state(F,:ntrec),2)] .= 0 + s3 = srcillum(J; time_mask) close(F) @test s1 ≈ s2 + @test norm(s3) < norm(s2) @test maximum(s1) > 0 end diff --git a/test/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl b/test/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl index b965ccc..572a5f7 100644 --- a/test/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl +++ b/test/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl @@ -184,8 +184,13 @@ end J = jacobian(F, m₀); s1 = srcillum(F, m₀); s2 = srcillum(J); + + time_mask = ones(Float32, state(F, :ntrec)) + time_mask[1:div(state(F,:ntrec),2)] .= 0 + s3 = srcillum(J; time_mask) close(F) @test s1 ≈ s2 + @test norm(s3) < norm(s2) @test maximum(s1) > 0 end diff --git a/test/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl b/test/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl index 2066c84..bcf1934 100644 --- a/test/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl +++ b/test/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl @@ -194,8 +194,13 @@ end J = jacobian(F, m₀); s1 = srcillum(F, m₀); s2 = srcillum(J); + + time_mask = ones(Float32, state(F, :ntrec)) + time_mask[1:div(state(F,:ntrec),2)] .= 0 + s3 = srcillum(J; time_mask) close(F) @test s1 ≈ s2 + @test norm(s3) < norm(s2) @test maximum(s1) > 0 end diff --git a/test/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl b/test/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl index 42789f0..64e3c67 100644 --- a/test/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl +++ b/test/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl @@ -188,8 +188,13 @@ end J = jacobian(F, m₀); s1 = srcillum(F, m₀); s2 = srcillum(J); + + time_mask = ones(Float32, state(F, :ntrec)) + time_mask[1:div(state(F,:ntrec),2)] .= 0 + s3 = srcillum(J; time_mask) close(F) @test s1 ≈ s2 + @test norm(s3) < norm(s2) @test maximum(s1) > 0 end diff --git a/test/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl b/test/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl index fd8f1af..9b9ea7f 100644 --- a/test/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl +++ b/test/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl @@ -186,8 +186,13 @@ end J = jacobian(F, m₀); s1 = srcillum(F, m₀); s2 = srcillum(J); + + time_mask = ones(Float32, state(F, :ntrec)) + time_mask[1:div(state(F,:ntrec),2)] .= 0 + s3 = srcillum(J; time_mask) close(F) @test s1 ≈ s2 + @test norm(s3) < norm(s2) @test maximum(s1) > 0 end