Open kerim371 opened 1 year ago
Oh, it seems negative spacing doesn't work:
ERROR: PyError ($(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))))) <class 'ValueError'>
ValueError('`nt` must be > 0')
File "/home/kerim/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/pysource/interface.py", line 46, in forward_rec
rec, _, I, _ = forward(model, src_coords, rec_coords, wavelet, save=False,
File "/home/kerim/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/pysource/propagators.py", line 30, in forward
src, rcv = src_rec(model, u, src_coords, rcv_coords, wavelet, nt)
File "/home/kerim/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/pysource/geom_utils.py", line 14, in src_rec
src = PointSource(name="src%s" % namef, grid=model.grid, ntime=nt,
File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/types/basic.py", line 873, in __new__
newobj._shape = cls.__shape_setup__(**kwargs)
File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/types/sparse.py", line 336, in __shape_setup__
raise ValueError('`nt` must be > 0')
Stacktrace:
[1] pyerr_check
@ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/exception.jl:62 [inlined]
[2] pyerr_check
@ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/exception.jl:66 [inlined]
[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}}, 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}}, kwargs::Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:fw, :space_order, :f0, :illum), Tuple{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, ::Vararg{Any, N} where N; kwargs::Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:fw, :space_order, :f0, :illum), Tuple{Bool, Int64, Float32, Bool}}})
@ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/pyfncall.jl:80
[11] (::JUDI.var"#189#190"{Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:fw, :space_order, :f0, :illum), Tuple{Bool, Int64, Float32, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, DataType})()
@ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/python_interface.jl:20
[12] (::JUDI.var"#1#2"{JUDI.var"#189#190"{Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:fw, :space_order, :f0, :illum), Tuple{Bool, Int64, Float32, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, 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"#189#190"{Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:fw, :space_order, :f0, :illum), Tuple{Bool, Int64, Float32, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, DataType}}, l::ReentrantLock)
@ Base ./lock.jl:187
[14] pylock(f::JUDI.var"#189#190"{Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:fw, :space_order, :f0, :illum), Tuple{Bool, Int64, Float32, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, DataType})
@ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:65
[15] wrapcall_data(::PyCall.PyObject, ::PyCall.PyObject, ::Vararg{Any, N} where N; kw::Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:fw, :space_order, :f0, :illum), Tuple{Bool, Int64, Float32, Bool}}})
@ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/python_interface.jl:19
[16] devito_interface(modelPy::PyCall.PyObject, srcGeometry::GeometryIC{Float32}, srcData::Matrix{Float32}, recGeometry::GeometryIC{Float32}, recData::Nothing, dm::Nothing, options::JUDIOptions, illum::Bool, fw::Bool)
@ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/python_interface.jl:72
[17] time_modeling(model_full::Model, srcGeometry::GeometryOOC{Float32}, srcData::Matrix{Float32}, recGeometry::GeometryOOC{Float32}, recData::Nothing, 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::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"{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::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] *
@ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/LinearOperators/operators.jl:173 [inlined]
[22] afoldl
@ ./operators.jl:533 [inlined]
[23] *(a::judiProjection{Float32}, b::judiModeling{Float32, :forward}, c::JUDI.jAdjoint{judiProjection{Float32}}, xs::judiVector{Float32, Matrix{Float32}})
@ Base ./operators.jl:560
[24] top-level scope
@ /mnt/HDD_2TB/130406/inversion/scripts/tests_modeling.jl:129
I'm actually not quite sure what value people usually put in the SegY header for these type of geometry, is there some open data you are using that shows what these values usually are?
Negative spacing will be very problematic yeah, negative origin would be the way to go as long as the segy is read correctly and put negative values for the rec depth.
@mloubout I'm not sure about open source data but in SEGY header positive elevations usually mean above sea level. Even though seismic software may treat this differently and I always look for the documentation whether above sea level is positive or negative.
With JUDI it seems the simplest solution is to treat RecGroupElevation
and SourceSurfaceElevation
without making hem absolute.
In this way if we have some registered data above sea level then we can set vertical origin to negative value, vertical spacing is positive and RecGroupElevation SourceSurfaceElevation
in SEGY header are negative.
I propose to add an argument to JUDI::Geometry
to multiply elevations.
For example make something like:
function Geometry(data::SegyIO.SeisBlock; key="source", segy_depth_key="", dt=nothing, t=nothing, elev_mult=1f0)
and use elev_mult
as multiplier to elevations.
It is worth to add information to the documentation that JUDI treats negative elevations above sea level and positive elevations below sea level to make it consistent with model origin/spacings.
I don't think adding this keyword to Geometry is a good idea, I tend to prefer to have a robust reader rather than having a lot of arguments. So if the abs
is the issue I would just remove it
Yes removing abs
would help.
But adding this information to the documentation is necessary as it is not easy to guess.
Hi,
Lets suppose we work with land data. Of course we have to take into account the relief (elevations). If the land elevation is above sea level (say
+100
meters) then probably model's vertical origin should start with-100
right? and vertical spacing is positive (say+25
m).The problem is that I don't know what sign should have
RecGroupElevation
andSourceSurfaceElevation
. I think they should be-100
m but anyway JUDI reads them by absolute value: link1, link2. In this case most likely the JUDI will treat them as they located below sea level.Also there is more difficult case when elevation on the profile vary from positive to negative values. Taking absolute elevations will lead to incorrect geometry.
If all my elevations are above sea level may I set vertical origin to positive (
+100
m) and vertical spacing to negative (-25
m)?