slimgroup / JUDI.jl

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

Dimension mismatch in `process_input_data` #211

Closed kerim371 closed 8 months ago

kerim371 commented 8 months ago

Hi,

I have a strange behaviour. When performing RTM with with OOC geometry on my local PC Ubuntu 20.04 there is no any problems. But when running it on the cloud CentOS 7 I get error:

ERROR: type GeometryOOC has no field xloc
Stacktrace:
  [1] getproperty(x::GeometryOOC{Float32}, f::Symbol)
    @ Base ./Base.jl:33
  [2] process_input_data(input::Vector{Float32}, geometry::GeometryOOC{Float32})
    @ JUDI ~/.julia/packages/JUDI/JEsVr/src/TimeModeling/Utils/auxiliaryFunctions.jl:630

But I can override this problem by determining the number of sources with nsrc = get_nsrc(geometry). That leads to the next error:

ERROR: DimensionMismatch("new dimensions (220, 301, 121) must be consistent with array size 7745100")
Stacktrace:
  [1] (::Base.var"#throw_dmrsa#234")(dims::Tuple{Int64, Int64, Int64}, len::Int64)
    @ Base ./reshapedarray.jl:41
  [2] reshape
    @ ./reshapedarray.jl:45 [inlined]
  [3] reshape
    @ ./reshapedarray.jl:116 [inlined]
  [4] process_input_data(input::Vector{Float32}, geometry::GeometryOOC{Float32})
    @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/Utils/auxiliaryFunctions.jl:631
  [5] process_input_data(J::judiJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}, q::Vector{Float32})
    @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/LinearOperators/operators.jl:195
  [6] multi_src_propagate(F::judiJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}, q::Vector{Float32})
    @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:79
  [7] *
    @ ~/.julia/dev/JUDI/src/TimeModeling/LinearOperators/operators.jl:176 [inlined]
  [8] (::JOLI.var"#173#177"{JOLI.joLinearOperator{Float32, Float32}, judiJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}})(v1::Vector{Float32})
    @ JOLI ~/.julia/packages/JOLI/S3LJy/src/joAbstractLinearOperator/base_functions.jl:139
  [9] *(A::JOLI.joLinearOperator{Float32, Float32}, v::Vector{Float32})
    @ JOLI ~/.julia/packages/JOLI/S3LJy/src/joAbstractLinearOperator/base_functions.jl:191
 [10] *(J::JOLI.joLinearOperator{Float32, Float32}, x::judiVector{Float32, Matrix{Float32}})
    @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/Types/abstract.jl:86
 [11] *(::JOLI.joLinearOperator{Float32, Float32}, ::judiJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}, ::judiVector{Float32, Matrix{Float32}})
    @ Base ./operators.jl:560
 [12] top-level scope
    @ ~/shared/MY_FOLDER/proc/rtm.jl:135

caused by: PyError ($(Expr(:escape, :(ccall(#= /home/kerim/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'pytools.prefork.ExecError'>
ExecError("error invoking '/home/kerim/shared_app/gcc/10.2/bin/gcc -O3 -g -march=native -fPIC -Wall -std=c99 -shared -fopenmp /tmp/devito-jitcache-uid1000/8b3bd02a06798a5b5f4db3d3ae57ef228af81907.c -lm -o /tmp/devito-jitcache-uid1000/8b3bd02a06798a5b5f4db3d3ae57ef228af81907.so': [Errno 12] Cannot allocate memory")
  File "/home/kerim/.julia/dev/JUDI/src/pysource/interface.py", line 359, in J_adjoint
    return J_adjoint_standard(model, src_coords, wavelet, rec_coords, recin,
  File "/home/kerim/.julia/dev/JUDI/src/pysource/interface.py", line 486, in J_adjoint_standard
    g, Iv, _ = gradient(model, residual, rec_coords, u, space_order=space_order, ic=ic,
  File "/home/kerim/.julia/dev/JUDI/src/pysource/propagators.py", line 129, in gradient
    summary = op(**kw)
  File "/home/kerim/.local/lib/python3.9/site-packages/devito/operator/operator.py", line 753, in __call__
    return self.apply(**kwargs)
  File "/home/kerim/.local/lib/python3.9/site-packages/devito/operator/operator.py", line 824, in apply
    cfunction = self.cfunction
  File "/home/kerim/.local/lib/python3.9/site-packages/devito/operator/operator.py", line 706, in cfunction
    self._jit_compile()
  File "/home/kerim/.local/lib/python3.9/site-packages/devito/operator/operator.py", line 691, in _jit_compile
    recompiled, src_file = self._compiler.jit_compile(self._soname,
  File "/home/kerim/.local/lib/python3.9/site-packages/devito/arch/compiler.py", line 360, in jit_compile
    _, _, _, recompiled = compile_from_string(self, target, code, src_file,
  File "/home/kerim/.local/lib/python3.9/site-packages/codepy/jit.py", line 439, in compile_from_string
    toolchain.build_extension(ext_file, source_paths, debug=debug)
  File "/home/kerim/.local/lib/python3.9/site-packages/codepy/toolchain.py", line 205, in build_extension
    result = call(cc_cmdline)
  File "/home/kerim/.local/lib/python3.9/site-packages/pytools/prefork.py", line 213, in call
    return forker.call(cmdline, cwd)
  File "/home/kerim/.local/lib/python3.9/site-packages/pytools/prefork.py", line 26, in call
    raise ExecError("error invoking '{}': {}".format(" ".join(cmdline), e))

Stacktrace:
  [1] pyerr_check
    @ ~/.julia/packages/PyCall/ilqDX/src/exception.jl:75 [inlined]
  [2] pyerr_check
    @ ~/.julia/packages/PyCall/ilqDX/src/exception.jl:79 [inlined]
  [3] _handle_error(msg::String)
    @ PyCall ~/.julia/packages/PyCall/ilqDX/src/exception.jl:96
  [4] macro expansion
    @ ~/.julia/packages/PyCall/ilqDX/src/exception.jl:110 [inlined]
  [5] #107
    @ ~/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:43 [inlined]
  [6] disable_sigint
    @ ./c.jl:458 [inlined]
  [7] __pycall!
    @ ~/.julia/packages/PyCall/ilqDX/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 ~/.julia/packages/PyCall/ilqDX/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, Nothing, String, Bool, Int64, Float32, Bool}}})
    @ PyCall ~/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:11
 [10] #pycall#112
    @ ~/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:80 [inlined]
 [11] #4
    @ ~/.julia/dev/JUDI/src/JUDI.jl:82 [inlined]
 [12] (::JUDI.var"#1#2"{JUDI.var"#4#5"{Tuple{PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, 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, Nothing, String, Bool, Int64, Float32, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}}})()
    @ JUDI ~/.julia/dev/JUDI/src/JUDI.jl:74
 [13] lock(f::JUDI.var"#1#2"{JUDI.var"#4#5"{Tuple{PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, 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, Nothing, String, Bool, Int64, Float32, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}}}, l::ReentrantLock)
    @ Base ./lock.jl:187
 [14] pylock(f::JUDI.var"#4#5"{Tuple{PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, 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, Nothing, String, Bool, Int64, Float32, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}})
    @ JUDI ~/.julia/dev/JUDI/src/JUDI.jl:71
 [15] rlock_pycall(::PyCall.PyObject, ::Type{Tuple{PyCall.PyArray, PyCall.PyObject, 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, Nothing, String, Bool, Int64, Float32, Bool}}})
    @ JUDI ~/.julia/dev/JUDI/src/JUDI.jl:81
 [16] 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, Nothing, String, Bool, Int64, Float32, Bool}}})
    @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/Modeling/python_interface.jl:46
 [17] 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 ~/.julia/dev/JUDI/src/TimeModeling/Modeling/python_interface.jl:136
 [18] macro expansion
    @ ~/.julia/dev/JUDI/src/TimeModeling/Modeling/time_modeling_serial.jl:46 [inlined]
 [19] macro expansion
    @ ./timing.jl:287 [inlined]
 [20] macro expansion
    @ ~/.julia/dev/JUDI/src/JUDI.jl:141 [inlined]
 [21] time_modeling(model_full::JUDI.IsoModel{Float32, 2}, srcGeometry::GeometryOOC{Float32}, srcData::Matrix{Float32}, recGeometry::GeometryOOC{Float32}, recData::Matrix{Float32}, dm::Nothing, op::Symbol, options::JUDIOptions, fw::Bool)
    @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/Modeling/time_modeling_serial.jl:45
 [22] propagate(F::judiJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}, q::judiVector{Float32, Matrix{Float32}})
    @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:9
 [23] macro expansion
    @ ~/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:36 [inlined]
 [24] macro expansion
    @ ./timing.jl:287 [inlined]
 [25] macro expansion
    @ ~/.julia/dev/JUDI/src/JUDI.jl:141 [inlined]
 [26] run_and_reduce(func::Function, #unused#::Nothing, nsrc::Int64, arg_func::JUDI.var"#267#268"{judiJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}, judiVector{Float32, Matrix{Float32}}})
    @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:35
 [27] multi_src_propagate(F::judiJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}, q::judiVector{Float32, Matrix{Float32}})
    @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:85
 [28] *
    @ ~/.julia/dev/JUDI/src/TimeModeling/LinearOperators/operators.jl:173 [inlined]
 [29] (::JOLI.var"#173#177"{JOLI.joLinearOperator{Float32, Float32}, judiJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}})(v1::judiVector{Float32, Matrix{Float32}})
    @ JOLI ~/.julia/packages/JOLI/S3LJy/src/joAbstractLinearOperator/base_functions.jl:139
 [30] *(J::JOLI.joLinearOperator{Float32, Float32}, x::judiVector{Float32, Matrix{Float32}})
    @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/Types/abstract.jl:86
 [31] *(::JOLI.joLinearOperator{Float32, Float32}, ::judiJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}, ::judiVector{Float32, Matrix{Float32}})
    @ Base ./operators.jl:560
 [32] top-level scope
    @ ~/shared/MY_FOLDER/proc/rtm.jl:135

As I said when I run the same code on my local Ubuntu 20.04 with the same JUDI and Devito versions I don't have any troubles and moreover it seems that process_input_data isn't invoked at all.

JUDI 3.3.8 Devito 4.8.2

kerim371 commented 8 months ago

And the script is:

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

# Load starting model (mlog - slowness built with Vs from logs; mvsp - built from VSP)
fid = h5open(model_file, "r")
n, d, o = read(fid, "n", "d", "o")
m0 = read(fid, "m")
close(fid)

n = Tuple(Int64(i) for i in n)
d = Tuple(Float32(i) for i in d)
o = Tuple(Float32(i) for i in o)
rho0 = rho_from_slowness(m0)
model0 = Model(n, d, o, m0, rho=rho0, nb=100)

x = (o[1]:d[1]:o[1]+(n[1]-1)*d[1])./1000f0
z = (o[2]:d[2]:o[2]+(n[2]-1)*d[2])./1000f0

global seabed_ind = Int.(round.(seabed./d[2]))
model0.m[:,1:seabed_ind] .= (1/vwater)^2

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

timeD = src_geometry.t[1]
dtD = src_geometry.dt[1]

# setup wavelet
f0 = 0.01f0     # kHz
wavelet = ricker_wavelet(timeD, dtD, f0)
q = judiVector(src_geometry, wavelet)

############################################## RTM #################################################

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

# Right-hand preconditioners (model topmute)
idx_wb = find_water_bottom(reshape(model0.m, model0.n))
Tm = judiTopmute(model0.n, idx_wb, 10)  # Mute water column
S = judiDepthScaling(model0)
Mr = S*Tm

# Left-hand side preconditioners
Ml = judiDataMute(q.geometry, d_obs.geometry, vp=1300f0, t0=0.1f0, mode=:reflection) # keep reflections

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

shot_from = 1
shot_to = length(d_obs)
shot_step = 1   # only compute each nth shot (that is enough)

indsrc = rand(shot_from:shot_from+shot_step-1):shot_step:shot_to

# Topmute
d_obs = Ml[indsrc]*d_obs[indsrc]
kerim371 commented 8 months ago

Oh I'm sorry I just figured out that the problem was insufficient RAM on the cloud machine. Once I increased the memory the problem was solved.