add tests for t0

src/TimeModeling/TimeModeling.jl

@@ -29,6 +29,7 @@ const dmType{T} = Union{Vector{T}, PhysicalParameter{T}}
 #############################################################################
 # Utils
 include("Utils/auxiliaryFunctions.jl")
+include("Utils/time_utilities.jl")
 include("Modeling/losses.jl")
 
 #############################################################################

src/TimeModeling/Types/GeometryStructure.jl

@@ -182,6 +182,8 @@ end
 
 Geometry(xloc::CoordT, yloc::CoordT, zloc::CoordT, dt::Vector{T}, nt::Vector{<:Integer}, t::Vector{T}) where {T<:Real} = GeometryIC{T}(xloc,yloc,zloc,dt,nt, t)
 Geometry(xloc::CoordT, yloc::CoordT, zloc::CoordT, dt::Vector{T}, nt::Vector{T}, t::Vector{T}) where {T<:Real} = GeometryIC{T}(xloc,yloc,zloc,dt,convert(Vector{Int64}, nt), t)
+Geometry(xloc::CoordT, yloc::CoordT, zloc::CoordT, t::StepRangeLen{T}) where {T<:Real} = GeometryIC{T}(xloc,yloc,zloc,[t])
+Geometry(xloc::CoordT, yloc::CoordT, zloc::CoordT, t::Vector{<:StepRangeLen{T}}) where {T<:Real} = GeometryIC{T}(xloc,yloc,zloc,t)
 
 # For easy 2D setup
 Geometry(xloc, zloc; kw...) = Geometry(xloc, 0 .* xloc, zloc; kw...)
@@ -236,7 +238,7 @@ function timings_from_segy(data, dt=nothing, t=nothing, t0=nothing)
 
     if isnothing(t)
         ntCell = _get_p(nothing, data, nsrc, "ns", Int, 1)
-        tCell = Float32.((ntCell .- 1) .* dtCell)
+        tCell = Float32.((ntCell .- 1) .* dtCell .+ t0)
     else
         tCell = as_src_list(t, nsrc)
     end

src/TimeModeling/Types/judiVector.jl

@@ -123,9 +123,9 @@ function judiVector(geometry::Geometry, data::SegyIO.SeisCon)
     return judiVector{Float32, SegyIO.SeisCon}(nsrc, geometry,dataCell)
 end
 
-judiVector(data::SegyIO.SeisBlock; segy_depth_key="RecGroupElevation", t=nothing) = judiVector(Geometry(data; key="receiver", segy_depth_key=segy_depth_key, t=t), data)
-judiVector(data::SegyIO.SeisCon; segy_depth_key="RecGroupElevation", t=nothing)= judiVector(Geometry(data; key="receiver", segy_depth_key=segy_depth_key, t=t), data)
-judiVector(data::Vector{SegyIO.SeisCon}; segy_depth_key="RecGroupElevation", t=nothing) = judiVector(Geometry(data; key="receiver", segy_depth_key=segy_depth_key, t=t), data)
+judiVector(data::SegyIO.SeisBlock; kw...) = judiVector(Geometry(data; key="receiver", kw...), data)
+judiVector(data::SegyIO.SeisCon; kw...)= judiVector(Geometry(data; key="receiver", kw...), data)
+judiVector(data::Vector{SegyIO.SeisCon}; kw...) = judiVector(Geometry(data; key="receiver", kw...), data)
 judiVector(geometry::Geometry, data::Vector{SegyIO.SeisCon}) =  judiVector{Float32, SegyIO.SeisCon}(length(data), geometry, data)
 
 ############################################################

src/TimeModeling/Utils/auxiliaryFunctions.jl

@@ -7,7 +7,7 @@
 
 export ricker_wavelet, get_computational_nt, calculate_dt, setup_grid, setup_3D_grid
 export convertToCell, limit_model_to_receiver_area, remove_out_of_bounds_receivers, extend_gradient
-export time_resample, remove_padding, subsample, process_input_data
+export remove_padding, subsample, process_input_data
 export generate_distribution, select_frequencies
 export devito_model, pad_sizes, pad_array
 export transducer, as_vec
@@ -490,80 +490,6 @@ function setup_3D_grid(xrec::Vector{Any}, yrec::Vector{Any}, zrec::Vector{Any})
     setup_3D_grid(xrec, yrec, zrec)
 end
 
-"""
-    time_resample(data, geometry_in, dt_new)
-
-Resample the input data with sinc interpolation from the current time sampling (geometrty_in) to the
-new time sampling `dt_new`.
-
-Parameters
-* `data`: Data to be reampled. If data is a matrix, resamples each column.
-* `geometry_in`: Geometry on which `data` is defined.
-* `dt_new`: New time sampling rate to interpolate onto.
-"""
-function time_resample(data::AbstractArray{T, N}, G_in::Geometry, dt_new::Real) where {T<:Real, N}
-    tend = step(G_in.taxis[1])*(size(data, 1) - 1) + first(G_in.taxis[1])
-    new_t = first(G_in.taxis[1]):dt_new:tend
-    return time_resample(data, G_in.taxis[1], new_t)
-end
-
-"""
-    time_resample(data, dt_in, dt_new)
-
-Resample the input data with sinc interpolation from the current time sampling dt_in to the 
-new time sampling `dt_new`.
-
-Parameters
-* `data`: Data to be reampled. If data is a matrix, resamples each column.
-* `dt_in`: Time sampling of input
-* `dt_new`: New time sampling rate to interpolate onto.
-"""
-function time_resample(data::AbstractArray{T, N}, t_in::StepRangeLen, t_new::StepRangeLen) where {T<:Real, N}
-    dt_in, dt_new = step(t_in), step(t_new)
-    if dt_new==dt_in
-        return data
-    elseif (dt_new % dt_in) == 0
-        rate = Int64(div(dt_new, dt_in))
-        return _time_resample(data, rate)
-    else
-        @juditime "Data time sinc-interpolation" begin
-            dataInterp = Float32.(SincInterpolation(data, t_in, t_new))
-        end
-        return dataInterp
-    end
-end
-
-time_resample(data::AbstractArray{T, N}, dt_in::Number, dt_new::Number, t::Number) where {T<:Real, N} =
-    time_resample(data, 0:dt_in:(dt_in*ceil(t/dt_in)), 0:dt_new:(dt_new*ceil(t/dt_new)))
-
-
-"""
-    time_resample(data, dt_in, geometry_in)
-
-Resample the input data with sinc interpolation from the current time sampling (dt_in) to the
-new time sampling `geometry_out`.
-
-Parameters
-* `data`: Data to be reampled. If data is a matrix, resamples each column.
-* `geometry_out`: Geometry on which `data` is to be interpolated.
-* `dt_in`: Time sampling rate of the `data.`
-"""
-function time_resample(data::AbstractArray{T, N}, dt_in::Real, G_out::Geometry{T}) where {T<:Real, N}
-    currt = range(0f0, step=dt_in, length=size(data, 1))
-    return time_resample(data, currt, G_out.taxis[1])
-end
-
-function time_resample(data::AbstractArray{T, N}, dt_in::Real, G_in::Geometry{T}, G_out::Geometry{T}) where {T<:Real, N}
-    currt = range(first(G_in.taxis[1]), step=dt_in, length=size(data, 1))
-    return time_resample(data, currt, G_out.taxis[1])
-end
-
-_time_resample(data::Matrix{T}, rate::Integer) where T = data[1:rate:end, :]
-_time_resample(data::PermutedDimsArray{T, 2, (2, 1), (2, 1), Matrix{T}}, rate::Integer) where {T<:Real} = data.parent[:, 1:rate:end]'
-
-SincInterpolation(Y::Matrix{T}, S::StepRangeLen{T}, Up::StepRangeLen{T}) where T<:Real = sinc.( (Up .- S') ./ (S[2] - S[1]) ) * Y
-SincInterpolation(Y::PermutedDimsArray{T, 2, (2, 1), (2, 1), Matrix{T}}, S::StepRangeLen{T}, Up::StepRangeLen{T}) where T<:Real = (Y.parent * sinc.( (Up' .- S) ./ (S[2] - S[1]) ))'
-
 """
     generate_distribution(x; src_no=1)
 
@@ -806,26 +732,3 @@ function filter_none(args::Tuple)
 end
 
 filter_none(x) = x
-
-
-"""
-    _maybe_pad_t0(q, qGeom, data, dataGeom)
-
-Pad zeros for data with non-zero t0, usually from a segy file so that time axis and array size match for the source and data.
-"""
-function _maybe_pad_t0(qIn::Matrix{T}, qGeom::Geometry, dObserved::Matrix{T}, dataGeom::Geometry) where T<:Number
-    if size(dObserved, 1) != size(qIn, 1)
-        dsize = size(qIn, 1) - size(dObserved, 1)
-        dt0 = get_t0(dataGeom, 1) - get_t0(qGeom, 1)
-        @assert dt0 != 0 && sign(dsize) == sign(dt0)
-        if dt0 > 0
-            dObserved = vcat(zeros(T, dsize, size(dObserved, 2)), dObserved)
-        else
-            qIn = vcat(zeros(T, -dsize,  size(qIn, 2)), qIn)
-        end
-    end
-    return qIn, dObserved
-end
-
-_maybe_pad_t0(qIn::judiVector{T, AT}, dObserved::judiVector{T, AT}) where{T<:Number, AT} =
-    _maybe_pad_t0(qIn.data, qIn.geometry, dObserved.data, dObserved.geometry)

src/TimeModeling/Utils/time_utilities.jl

@@ -0,0 +1,140 @@
+export time_resample
+
+
+"""
+    time_resample(data, geometry_in, dt_new)
+
+Resample the input data with sinc interpolation from the current time sampling (geometrty_in) to the
+new time sampling `dt_new`.
+
+Parameters
+* `data`: Data to be reampled. If data is a matrix, resamples each column.
+* `geometry_in`: Geometry on which `data` is defined.
+* `dt_new`: New time sampling rate to interpolate onto.
+"""
+function time_resample(data::AbstractArray{T, N}, G_in::Geometry, dt_new::Real) where {T<:Real, N}
+    tend = step(G_in.taxis[1])*(size(data, 1) - 1) + first(G_in.taxis[1])
+    new_t = first(G_in.taxis[1]):dt_new:tend
+    return time_resample(data, G_in.taxis[1], new_t)
+end
+
+"""
+    time_resample(data, dt_in, dt_new)
+
+Resample the input data with sinc interpolation from the current time sampling dt_in to the 
+new time sampling `dt_new`.
+
+Parameters
+* `data`: Data to be reampled. If data is a matrix, resamples each column.
+* `dt_in`: Time sampling of input
+* `dt_new`: New time sampling rate to interpolate onto.
+"""
+function time_resample(data::AbstractArray{T, N}, t_in::StepRangeLen, t_new::StepRangeLen) where {T<:Real, N}
+    dt_in, dt_new = step(t_in), step(t_new)
+    if dt_new==dt_in
+        return data
+    elseif (dt_new % dt_in) == 0
+        rate = Int64(div(dt_new, dt_in))
+        return _time_resample(data, rate)
+    else
+        @juditime "Data time sinc-interpolation" begin
+            dataInterp = Float32.(SincInterpolation(data, t_in, t_new))
+        end
+        return dataInterp
+    end
+end
+
+time_resample(data::AbstractArray{T, N}, dt_in::Number, dt_new::Number, t::Number) where {T<:Real, N} =
+    time_resample(data, 0:dt_in:(dt_in*ceil(t/dt_in)), 0:dt_new:(dt_new*ceil(t/dt_new)))
+
+
+"""
+    time_resample(data, dt_in, geometry_in)
+
+Resample the input data with sinc interpolation from the current time sampling (dt_in) to the
+new time sampling `geometry_out`.
+
+Parameters
+* `data`: Data to be reampled. If data is a matrix, resamples each column.
+* `geometry_out`: Geometry on which `data` is to be interpolated.
+* `dt_in`: Time sampling rate of the `data.`
+"""
+function time_resample(data::AbstractArray{T, N}, dt_in::Real, G_out::Geometry{T}) where {T<:Real, N}
+    currt = range(0f0, step=dt_in, length=size(data, 1))
+    return time_resample(data, currt, G_out.taxis[1])
+end
+
+function time_resample(data::AbstractArray{T, N}, dt_in::Real, G_in::Geometry{T}, G_out::Geometry{T}) where {T<:Real, N}
+    t0 = min(get_t0(G_in, 1), get_t0(G_out, 1))
+    currt = range(t0, step=dt_in, length=size(data, 1))
+    return time_resample(data, currt, G_out.taxis[1])
+end
+
+_time_resample(data::Matrix{T}, rate::Integer) where T = data[1:rate:end, :]
+_time_resample(data::PermutedDimsArray{T, 2, (2, 1), (2, 1), Matrix{T}}, rate::Integer) where {T<:Real} = data.parent[:, 1:rate:end]'
+
+SincInterpolation(Y::Matrix{T}, S::StepRangeLen{T}, Up::StepRangeLen{T}) where T<:Real = sinc.( (Up .- S') ./ (S[2] - S[1]) ) * Y
+SincInterpolation(Y::PermutedDimsArray{T, 2, (2, 1), (2, 1), Matrix{T}}, S::StepRangeLen{T}, Up::StepRangeLen{T}) where T<:Real = (Y.parent * sinc.( (Up' .- S) ./ (S[2] - S[1]) ))'
+
+
+"""
+    _maybe_pad_t0(q, qGeom, data, dataGeom)
+
+Pad zeros for data with non-zero t0, usually from a segy file so that time axis and array size match for the source and data.
+"""
+function _maybe_pad_t0(qIn::Matrix{T}, qGeom::Geometry, dObserved::Matrix{T}, dataGeom::Geometry) where T<:Number
+    dt0 = get_t0(dataGeom, 1) - get_t0(qGeom, 1)
+    Dt = get_t(dataGeom, 1) - get_t(qGeom, 1)
+    dsize = size(qIn, 1) - size(dObserved, 1)
+    # Same times, do nothing
+    if dsize == 0 && dt0 == 0 && Dt == 0
+        return qIn, dObserved
+    # First case, same size, then it's a shift
+    elseif dsize == 0 && dt0 != 0 && Dt != 0
+        # Shift means both t0 and t same sign difference
+        @assert sign(dt0) == sign(Dt)
+        pad_size = Int(div(get_t0(dataGeom, 1), get_dt(dataGeom, 1)))
+        if dt0 > 0
+            # Data has larger t0, pad data left and q right
+            dObserved = vcat(zeros(T, pad_size, size(dObserved, 2)), dObserved)
+            qIn = vcat(qIn, zeros(T, pad_size, size(qIn, 2)))
+        else
+            # q has larger t0, pad data right and q left
+            dObserved = vcat(dObserved, zeros(T, pad_size, size(dObserved, 2)))
+            qIn = vcat(zeros(T, pad_size, size(qIn, 2)), qIn)
+        end
+    elseif dsize !=0
+        # We might still have differnt t0 and t
+        # Pad so that we go from smaller dt to largest t
+        ts = min(get_t0(qGeom, 1), get_t0(dataGeom, 1))
+        te = max(get_t(qGeom, 1), get_t(dataGeom, 1))
+
+        pdatal = Int(div(get_t0(dataGeom, 1) - ts, get_dt(dataGeom, 1)))
+        pdatar = Int(div(te - get_t(dataGeom, 1), get_dt(dataGeom, 1)))
+        dObserved = vcat(zeros(T, pdatal, size(dObserved, 2)), dObserved, zeros(T, pdatar, size(dObserved, 2)))
+
+        pql = Int(div(get_t0(qGeom, 1) - ts, get_dt(qGeom, 1)))
+        pqr = Int(div(te - get_t(qGeom, 1), get_dt(qGeom, 1)))
+        qIn = vcat(zeros(T, pql, size(qIn, 2)), qIn, zeros(T, pqr, size(qIn, 2)))
+    else
+        throw(judiMultiSourceException("""
+            Data and source have different
+            t0 : $((get_t0(dataGeom, 1), get_t0(qGeom, 1)))
+            and t: $((get_t(dataGeom, 1), get_t(qGeom, 1)))
+            and are not compatible in size for padding: $((size(qIn, 1), size(dObserved, 1)))"""))
+    end
+    return qIn, dObserved
+end
+
+pad_msg = """
+    This is an internal method for single source propatation,
+    only single-source judiVectors are supported
+"""
+
+function _maybe_pad_t0(qIn::judiVector{T, Matrix{T}}, dObserved::judiVector{T, Matrix{T}}) where{T<:Number}
+    @assert qIn.nsrc == 1 || throw(judiMultiSourceException(pad_msg))
+    return _maybe_pad_t0(qIn.data[1], qIn.geometry[1], dObserved.data[1], dObserved.geometry[1])
+end
+
+_maybe_pad_t0(qIn::judiVector{T, AT}, dObserved::judiVector{T, AT}) where{T<:Number, AT} =
+    _maybe_pad_t0(get_data(qIn), get_data(dObserved))

src/pysource/fields_exprs.py

@@ -76,6 +76,7 @@ def extented_src(model, weight, wavelet, q=0):
         return q
     ws, wt = lr_src_fields(model, weight, wavelet)
     if model.is_tti:
+        q = (0, 0) if q == 0 else q
         return (q[0] + ws * wt, q[1] + ws * wt)
     return q + ws * wt
 

src/pysource/interface.py

@@ -553,7 +553,7 @@ def J_adjoint_checkpointing(model, src_coords, wavelet, rec_coords, recin, space
 
     nt = wavelet.shape[0]
     rec = Receiver(name='rec', grid=model.grid, ntime=nt, coordinates=rec_coords)
-    kwg['srcv'] = rec
+    kwg['srcv1' if model.is_tti else 'srcv'] = rec
     # Wavefields to checkpoint
     cpwf = [uu for uu in as_tuple(u)]
     if model.is_viscoacoustic:

src/pysource/kernels.py

@@ -164,9 +164,9 @@ def tti_kernel(model, u1, u2, fw=True, q=None):
 
     if 'nofsdomain' in model.grid.subdomains:
         # Water at free surface, no anisotrpy
-        acout_ttip = [Eq(u1_n, stencilp.subs(model.zero_thomsen))]
-        acout_ttir = [Eq(u2_n, stencilr.subs(model.zero_thomsen))]
-        pdea = freesurface(model, (acout_ttip, acout_ttir), (u1, u2))
+        acout_ttip = Eq(u1_n, stencilp.subs(model.zero_thomsen))
+        acout_ttir = Eq(u2_n, stencilr.subs(model.zero_thomsen))
+        pdea = freesurface(model, [acout_ttip, acout_ttir], (u1, u2))
         # Standard PDE in subsurface
         first_stencil = Eq(u1_n, stencilp, subdomain=model.grid.subdomains['nofsdomain'])
         second_stencil = Eq(u2_n, stencilr, subdomain=model.grid.subdomains['nofsdomain'])