search

slimgroup / JUDI.jl

Julia Devito inversion.
https://slimgroup.github.io/JUDI.jl
MIT License
96 stars 29 forks source link

RTM in frequency domain result in error #164

Closed kerim371 closed 1 year ago

kerim371 commented 1 year ago

Hi,

I tried to perform RTM in frequency domain but I got an error:

container = segy_scan(prestk_dir, prestk_file, ["SourceX", "SourceY", "GroupX", "GroupY", "RecGroupElevation", "SourceSurfaceElevation", "dt"])
d_obs = judiVector(container; segy_depth_key = "RecGroupElevation")

# Set up wavelet and source vector
src_geometry = Geometry(container; key = "source", segy_depth_key = "SourceSurfaceElevation")
q = judiVector(src_geometry, wavelet)

t_sub = 2

# JUDI options
jopt = JUDI.Options(
    space_order=16,
    limit_m = true,
    buffer_size = buffer_size,
    optimal_checkpointing=false,
    IC = "isic")

# Setup operators
Pr = judiProjection(d_obs.geometry)
F = judiModeling(model0; options=jopt)
Ps = judiProjection(q.geometry)
J = judiJacobian(Pr*F*adjoint(Ps), q)

q_dist = generate_distribution(q)
# Set up random frequencies
nfreq = 20
J.options.frequencies = Array{Any}(undef, d_obs.nsrc)
J.options.dft_subsampling_factor = 8
for k=1:d_obs.nsrc
    J.options.frequencies[k] = select_frequencies(q_dist; fmin=0.003, fmax=0.07, nf=nfreq)
end

# Add random noise
d_lin = get_data(d_obs)
for j=1:d_lin.nsrc
    noise = randn(Float32, size(d_lin.data[j]))
    d_lin.data[j] += (noise/norm(vec(noise), Inf)*0.02f0)
end

# Topmute
d_lin = Ml_ref*d_lin

# RTM
rtm = adjoint(J)*d_lin

and the error:

/tmp/devito-jitcache-uid1000/86c2567d795373c88a9e10272fe3adb595c1c133.c: In function ‘gradient’:
/tmp/devito-jitcache-uid1000/86c2567d795373c88a9e10272fe3adb595c1c133.c:164:28: error: ‘r4’ undeclared (first use in this function); did you mean ‘r14’?
  164 |               float r14 = -r4*ufru[freq_dim][x + 1][y + 1];
      |                            ^~
      |                            r14
/tmp/devito-jitcache-uid1000/86c2567d795373c88a9e10272fe3adb595c1c133.c:164:28: note: each undeclared identifier is reported only once for each function it appears in
/tmp/devito-jitcache-uid1000/86c2567d795373c88a9e10272fe3adb595c1c133.c:165:27: error: ‘r5’ undeclared (first use in this function); did you mean ‘r15’?
  165 |               float r13 = r5*ufiu[freq_dim][x + 1][y + 1];
      |                           ^~
      |                           r15
FAILED compiler invocation:gcc-11 -O3 -g -fPIC -Wall -std=c99 -march=native -Wno-unused-result -Wno-unused-variable -Wno-unused-but-set-variable -ffast-math -shared -fopenmp /tmp/devito-jitcache-uid1000/86c2567d795373c88a9e10272fe3adb595c1c133.c -lm -o /tmp/devito-jitcache-uid1000/86c2567d795373c88a9e10272fe3adb595c1c133.so
/tmp/devito-jitcache-uid1000/3fc5d91317078abb08f733303d1646cbd30c2698.c: In function ‘gradient’:
/tmp/devito-jitcache-uid1000/3fc5d91317078abb08f733303d1646cbd30c2698.c:143:26: error: ‘r4’ undeclared (first use in this function); did you mean ‘r14’?
  143 |             float r14 = -r4*ufru[freq_dim][x + 1][y + 1];
      |                          ^~
      |                          r14
/tmp/devito-jitcache-uid1000/3fc5d91317078abb08f733303d1646cbd30c2698.c:143:26: note: each undeclared identifier is reported only once for each function it appears in
/tmp/devito-jitcache-uid1000/3fc5d91317078abb08f733303d1646cbd30c2698.c:144:25: error: ‘r5’ undeclared (first use in this function); did you mean ‘r15’?
  144 |             float r13 = r5*ufiu[freq_dim][x + 1][y + 1];
      |                         ^~
      |                         r15
FAILED compiler invocation:gcc-11 -O3 -g -fPIC -Wall -std=c99 -march=native -Wno-unused-result -Wno-unused-variable -Wno-unused-but-set-variable -ffast-math -shared -fopenmp /tmp/devito-jitcache-uid1000/3fc5d91317078abb08f733303d1646cbd30c2698.c -lm -o /tmp/devito-jitcache-uid1000/3fc5d91317078abb08f733303d1646cbd30c2698.so
ERROR: /PyErrorH ($(Expr(:escape, :(ccall(#= /home/kerim/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw)))))D D<class 'codepy.CompileError'>_
2CompileError('module compilation failed')T
B  File "/home/kerim/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/pysource/interface.py", line 353, in J_adjoint
    return J_adjoint_freq(model, src_coords, wavelet, rec_coords, recin, ws=ws,
/  File "/home/kerim/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/pysource/interface.py", line 424, in J_adjoint_freq
    g, Iv, _ = gradient(model, residual, rec_coords, u, space_order=space_order, ic=ic,
V  File "/home/kerim/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/pysource/propagators.py", line 127, in gradient
    summary = op(**kw)
i  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/operator/operator.py", line 686, in __call__
    return self.apply(**kwargs)
k  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/operator/operator.py", line 757, in apply
    cfunction = self.cfunction
i  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/operator/operator.py", line 639, in cfunction
    self._jit_compile()
n  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/operator/operator.py", line 624, in _jit_compile
    recompiled, src_file = self._compiler.jit_compile(self._soname,
g  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/arch/compiler.py", line 319, in jit_compile
    _, _, _, recompiled = compile_from_string(self, target, code, src_file,
G  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/codepy/jit.py", line 433, in compile_from_string
    toolchain.build_extension(ext_file, source_paths, debug=debug)
r  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/codepy/toolchain.py", line 210, in build_extension
    raise CompileError("module compilation failed")
a
Stacktrace:b
e n [1]/ lpyerr_checki
n    @ e~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/12exception.jl/:v62en [inlined]v
/ b [2]i npyerr_check/
julia> urce /mnt/HDD_2TB/VikingGraben/line12/venv/bin/activate

  [3] _handle_error(msg::String)
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/exception.jl:83
  [4] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/exception.jl:97 [inlined]
  [5] #107
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/pyfncall.jl:43 [inlined]
  [6] disable_sigint
    @ ./c.jl:458 [inlined]
  [7] __pycall!
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/pyfncall.jl:42 [inlined]
  [8] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, nargs::Int64, kw::PyCall.PyObject)
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/pyfncall.jl:29
  [9] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{10, Symbol}, NamedTuple{(:fw, :t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :dft_sub, :f0, :illum), Tuple{Bool, Int64, Int64, Bool, Vector{Float32}, String, Bool, Int64, Float32, Bool}}})
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/pyfncall.jl:11
 [10] pycall(::PyCall.PyObject, ::Type{Tuple{PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}}, ::PyCall.PyObject, ::Vararg{Any, N} where N; kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{10, Symbol}, NamedTuple{(:fw, :t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :dft_sub, :f0, :illum), Tuple{Bool, Int64, Int64, Bool, Vector{Float32}, String, Bool, Int64, Float32, Bool}}})
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/pyfncall.jl:80
 [11] (::JUDI.var"#198#199"{Base.Iterators.Pairs{Symbol, Any, NTuple{10, Symbol}, NamedTuple{(:fw, :t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :dft_sub, :f0, :illum), Tuple{Bool, Int64, Int64, Bool, Vector{Float32}, String, Bool, Int64, Float32, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, DataType})()
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/python_interface.jl:52
 [12] (::JUDI.var"#1#2"{JUDI.var"#198#199"{Base.Iterators.Pairs{Symbol, Any, NTuple{10, Symbol}, NamedTuple{(:fw, :t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :dft_sub, :f0, :illum), Tuple{Bool, Int64, Int64, Bool, Vector{Float32}, String, Bool, Int64, Float32, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, DataType}})()
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:68
 [13] lock(f::JUDI.var"#1#2"{JUDI.var"#198#199"{Base.Iterators.Pairs{Symbol, Any, NTuple{10, Symbol}, NamedTuple{(:fw, :t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :dft_sub, :f0, :illum), Tuple{Bool, Int64, Int64, Bool, Vector{Float32}, String, Bool, Int64, Float32, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, DataType}}, l::ReentrantLock)
    @ Base ./lock.jl:187
 [14] pylock(f::JUDI.var"#198#199"{Base.Iterators.Pairs{Symbol, Any, NTuple{10, Symbol}, NamedTuple{(:fw, :t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :dft_sub, :f0, :illum), Tuple{Bool, Int64, Int64, Bool, Vector{Float32}, String, Bool, Int64, Float32, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, DataType})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:65
 [15] wrapcall_grad(::PyCall.PyObject, ::PyCall.PyObject, ::Vararg{Any, N} where N; kw::Base.Iterators.Pairs{Symbol, Any, NTuple{10, Symbol}, NamedTuple{(:fw, :t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :dft_sub, :f0, :illum), Tuple{Bool, Int64, Int64, Bool, Vector{Float32}, String, Bool, Int64, Float32, Bool}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/python_interface.jl:51
 [16] devito_interface(modelPy::PyCall.PyObject, srcGeometry::GeometryIC{Float32}, srcData::Matrix{Float32}, recGeometry::GeometryIC{Float32}, recData::Matrix{Float32}, dm::Nothing, options::JUDIOptions, illum::Bool, fw::Bool)
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/python_interface.jl:143
 [17] time_modeling(model_full::Model, srcGeometry::GeometryOOC{Float32}, srcData::Matrix{Float32}, recGeometry::GeometryOOC{Float32}, recData::Matrix{Float32}, dm::Nothing, op::Symbol, options::JUDIOptions, fw::Bool)
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/time_modeling_serial.jl:39
 [18] propagate(F::judiJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}, q::judiVector{Float32, Matrix{Float32}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:9
 [19] run_and_reduce(func::Function, #unused#::Nothing, nsrc::Int64, arg_func::JUDI.var"#228#229"{judiJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}, judiVector{Float32, Matrix{Float32}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:32
 [20] multi_src_propagate(F::judiJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}, q::judiVector{Float32, Matrix{Float32}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:67
 [21] *(F::judiJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}, q::judiVector{Float32, Matrix{Float32}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/LinearOperators/operators.jl:173
 [22] top-level scope
    @ /mnt/HDD_2TB/VikingGraben/line12/inversion/scripts/rtm.jl:192
mloubout commented 1 year ago

That's quite odd this is tested. What version of devito do you use?

kerim371 commented 1 year ago

version 4.7.1

import devito
devito.__version__ 
'4.7.1'
mloubout commented 1 year ago

Looks like the issue is the combination of dft + dft_subsamp + isic will look into it, in the mean time turning off dft_subsamp should avoid the error

kerim371 commented 1 year ago

Right, without subsampling it works

mloubout commented 1 year ago

Was a bug in devito, will be fixed with the next devito release once https://github.com/devitocodes/devito/pull/2045 is merged

mloubout commented 1 year ago

If you switch to devito master it should work now, will be default with the next devito release.