adds a time mask for the src illumination integration

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

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

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

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

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

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

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

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

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
 

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
 

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
 

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
 

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
 

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