FWI: Dimension mismatch when limit_m=true

[kerim371]

Hi,

Starting from v3.1.9 I'm no longer able to run FWI with limit_m=true with my field data while fwi_example_2D.jl works fine even if model is limited.

The error I get:

ERROR: DimensionMismatch("new dimensions (477450,) must be consistent with array size 117000")
Stacktrace:
 [1] (::Base.var"#throw_dmrsa#234")(dims::Tuple{Int64}, len::Int64)
   @ Base ./reshapedarray.jl:41
 [2] reshape(a::Matrix{Float32}, dims::Tuple{Int64})
   @ Base ./reshapedarray.jl:45
 [3] materialize!(A::PhysicalParameter{Float32}, B::Base.Broadcast.Broadcasted{Base.Broadcast.ArrayStyle{PhysicalParameter}, Nothing, typeof(identity), Tuple{PhysicalParameter{Float32}}})
   @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Types/ModelStructure.jl:223
 [4] multi_src_fg!(G::PhysicalParameter{Float32}, model::Model, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}, dm::Nothing; options::JUDIOptions, kw::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:nlind, :lin), Tuple{Bool, Bool}}})
   @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:88
 [5] #multi_exp_fg!#186
   @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:155 [inlined]
 [6] fwi_objective!(G::PhysicalParameter{Float32}, model::Model, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}; options::JUDIOptions, kw::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:136
 [7] fwi_objective(model::Model, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}; options::JUDIOptions, kw::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:101
 [8] top-level scope
   @ /mnt/HDD_2TB/VikingGraben/line12/inversion/scripts/fwi_acoustic.jl:90

The simplified code I use:

# JUDI options
buffer_size = 0f0    # limit model (meters)

# Load data and create data vector
block = segy_read(prestk_dir * prestk_file)
# container = segy_scan(prestk_dir, prestk_file, ["SourceX", "SourceY", "GroupX", "GroupY", "RecGroupElevation", "SourceSurfaceElevation", "dt"])
d_obs = judiVector(block; segy_depth_key = "RecGroupElevation")

srcx = Float32.(get_header(block, "SourceX")[:,1])
grpx = Float32.(get_header(block, "GroupX")[:,1])
min_x = minimum([minimum(srcx), minimum(grpx)])-buffer_size
max_x = maximum([maximum(srcx), maximum(grpx)])+buffer_size

# Load starting model
n, d, o, m0 = read(h5open(model_file, "r"), "n", "d", "o", "m")

# Trim model
m0x = o[1]:d[1]:o[1]+(size(m0)[1]-1)*d[1]
x_ind = [i for i in eachindex(m0x) if min_x-d[1] <= m0x[i] <= max_x+d[1]]
m0 = m0[x_ind,:]
n = (length(x_ind), n[2])
d = (d[1], d[2])
o = (o[1]+(x_ind[1]-1)*d[1], o[2])
model0 = Model(n, d, o, m0)

# Set up wavelet and source vector
src_geometry = Geometry(block; key = "source", segy_depth_key = "SourceSurfaceElevation")
f0 = 0.006     # kHz
wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], f0)
q = judiVector(src_geometry, wavelet)

# JUDI options
jopt = JUDI.Options(
    limit_m = true,
    buffer_size = buffer_size,
    optimal_checkpointing=false)

# Select batch and calculate gradient (set to a single shot to accelerate crash experiment)
fval, gradient = fwi_objective(model0, q[10], d_obs[10], options=jopt)

Here:
o = (4968.5, 0.0)
n = (1061, 450)
d = (12.5, 12.5)
min_x = 4975.0f0 (from segy source and rec positions)
max_x = 18212.0f0 (from segy source and rec positions)

Each shot in SEGY contains 120 receivers

In the error the following numbers may be treated as:
477450 is equal to *(n...)
117000/n[2] = 260

[mloubout]

Thanks for raising this. I think I know where the issue is and will be fixing asap.

Basically, as of 3.1.9 each shot gradient is only the size of the aperture and the PhysicalParameter arithmetic takes care of summing different ones. But I didn't think of the case where the summer size wouldn't match the total model so I need to tweak it.

[mloubout]

Is fixed now