Exception when limit_m=true with custom container

[kerim371]

Hi,

I use JUDI 2.6.8 with custom container (based on hdf5 instead of segy)

I noticed that I get exceptions when limitting model (limit_m=true and buffer_size=100 for example):

JULIA: JOLI.joUtilsException("type mismatch: expected Float32 instead of Any in RDT from *(judiPDEfull,judiVector): / Proj*F*Proj' / JUDI.judiPDEfull{Float32, Float32} / Any")
Stacktrace:
 [1] jo_check_type_match(DT1::DataType, DT2::DataType, where::String)
   @ JOLI ~/.julia/packages/JOLI/g80Jj/src/joUtils.jl:329
 [2] *(A::JUDI.judiPDEfull{Float32, Float32}, v::JUDI.judiVector{Float32, Matrix{Float32}})
   @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/LinearOperators/judiPDEfull.jl:108
 [3] afoldl
   @ ./operators.jl:533 [inlined]
 [4] *
   @ ./operators.jl:560 [inlined]
 [5] H5ForwardModeling(; model::JUDI.Model, opt::JUDI.Options, q::JUDI.judiVector{Float32, Matrix{Float32}}, recGeometry::ColadaSeismic.H5GeometryOOC)
   @ ColadaSeismic ~/.julia/dev/ColadaSeismic/src/TimeModeling/h5judi_forward_modeling.jl:21
 [6] H5Modeling(; h5vel::PyObject, h5density::Nothing, h5epsilon::Nothing, h5delta::Nothing, h5tetha::Nothing, h5phi::Nothing, model_xkey::String, model_ykey::String, h5geom::PyObject, geom_src_pkey::String, geom_src_xkey::String, geom_src_ykey::String, geom_src_zkey::String, geom_rec_xkey::String, geom_rec_ykey::String, geom_rec_zkey::String, opt::JUDI.Options, src_frq::Float64, vmin::Float64, vmax::Float64, computation_type::String, fwi_niter::Int64, fwi_batchsize::Int64, twri_opt::JUDI.TWRIOptions)
   @ ColadaSeismic ~/.julia/dev/ColadaSeismic/src/TimeModeling/h5judi_modeling.jl:210
 [7] invokelatest(::Any; kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{24, Symbol}, NamedTuple{(:h5vel, :h5density, :h5epsilon, :h5delta, :h5tetha, :h5phi, :model_xkey, :model_ykey, :h5geom, :geom_src_pkey, :geom_src_xkey, :geom_src_ykey, :geom_src_zkey, :geom_rec_xkey, :geom_rec_ykey, :geom_rec_zkey, :opt, :src_frq, :vmin, :vmax, :computation_type, :fwi_niter, :fwi_batchsize, :twri_opt), Tuple{PyObject, Nothing, Nothing, Nothing, Nothing, Nothing, String, String, PyObject, String, String, String, String, String, String, String, JUDI.Options, Float64, Float64, Float64, String, Int64, Int64, JUDI.TWRIOptions}}})
   @ Base ./essentials.jl:710
 [8] _pyjlwrap_call(f::Function, args_::Ptr{PyCall.PyObject_struct}, kw_::Ptr{PyCall.PyObject_struct})
   @ PyCall ~/.julia/packages/PyCall/7a7w0/src/callback.jl:32
 [9] pyjlwrap_call(self_::Ptr{PyCall.PyObject_struct}, args_::Ptr{PyCall.PyObject_struct}, kw_::Ptr{PyCall.PyObject_struct})
   @ PyCall ~/.julia/packages/PyCall/7a7w0/src/callback.jl:44>

To get some information what happens there I added some @info to TimeModeling/LinearOperators/judiPDEfull.jl:

function *(A::judiPDEfull{ADDT,ARDT}, v::judiVector{vDT, AT}) where {ADDT, ARDT,vDT, AT}
    A.n == size(v,1) || throw(judiPDEfullException("Shape mismatch: A:$(size(A)), v: $(size(v))"))
    if compareGeometry(A.srcGeometry,v.geometry) == false && compareGeometry(A.recGeometry,v.geometry) == false
        throw(judiPDEfullException("geometry mismatch"))
    end
    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiPDEfull,judiVector):",A.name,typeof(A),vDT]," / "))
    V = A.fop(v)
    @info "eltype(V): $(eltype(V))\n"
    @info "typeof(V): $(typeof(V))\n"
    @info "typeof(ARDT): $(typeof(ARDT))\n"
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDEfull,judiVector):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

and the output is:

[ Info: eltype(V): Any
[ Info: typeof(V): Distributed.RemoteException
[ Info: typeof(ARDT): DataType

But if I run modeling_basic_2D (and _3D) example then everything is fine and the output:

[ Info: eltype(V): Float32
[ Info: typeof(V): judiVector{Float32, SeisCon}
[ Info: typeof(ARDT): DataType

Maybe you may give a hint where can I look to find a reason for that exception?
If I don't use model limitting then the calculation (limit_m=false) seem to run fine.

[mloubout]

This RemoteException usually means the python devito part errorer and pycall returned an error message instead of the result. Does it error with standard judiVector as well (segy and normal ones)?

[kerim371]

Does it error with standard judiVector as well (segy and normal ones)?

Do you mean have I tried to run the modeling with same data but with judiVector constructed with explicit GeometryIC and judiVector constructed with SeisCon?
If so I have not tried this with my data and probably it should be doable pretty simple with GeometryIC but it would take some time to export my data from HDF5 to segy (it is small few megabytes but some code is required)

But as I said above I tried to run JUDI examples 2D and 3D and both of them run perfectly.

I remembered that my model has some non-zero origin. May this be the reason? In an one-two hours (now JUDI runs on my PC doing test modeling) I will try your recommendations and try to set the model's origin to zero.

[mloubout]

I think it's fine without segy just the standard judi vector and see how it differs with yours. A non zero origin should not be an issue and should be tested as well

[kerim371]

@mloubout it seems I've found something.

Neither custom InCore Geometry nor Out Of Core Geometry works with limit_m=true and non-zero origin.
But if I set model origin (and shift geometry by the same value of course) to zero then both of them work.

In the same time JUDI 3D example works without any problem even if I set limit_m=true, buffer_size=100 and origin to any reasonable value.

To better explain the situation I'm testing Teapot Dome data.
It has some origin and it is rotated a little bit in XY plane.
To launch JUDI I calculate the angle and apply it to coordinates so that the model and geometry is XY alined. I do that at data loading time.
There should not be any problems with rotations as when I shifted the origin I still had to apply rotation. It should be fine.

[mloubout]

So looks like the limit_model_to_receiver_area function isn't working correctly with your geometry. You could go into it and check what the new size and origines are after crop to check all still make sense

[kerim371]

Hmm, I ecnountered another ptoblem.

I aborted computation and seems after that when I run modeling my app crashes as julia consumes all the RAM (about 12 Gb).
Don't you know why?

[mloubout]

Will depend on the size of your problem so may run out of memory yes.

[kerim371]

The model size 188 * 345 * 1501 (about 300-400 Mb),
The modeled SP seismogram (geometry) has 2500 traces and 1501 time sample.

Is it normal that the computation takes 13.3 Gb?
In the output I can see:

recompiling for non-existent cache dir (/tmp/devito-codepy-uid1000/5844cc5/55e55150988140bc4c14e71919a4c451).
recompiling for non-existent cache dir (/tmp/devito-codepy-uid1000/3ab2870/5df1815623306e5f3ecfa5bc29ad365c).

Also my PC has 8 cores (4 physical, 4 logical), Ubuntu 20.04.
I can see that computation takes exactly 50% of CPU.
Is there any reason why it can't be say 90% or is it somehow restricted? For example I remember by default WSL2 on Windows involes only half of computer power. Probably it is somehow similar? (P.S. I use Ubuntu, not WSL2, it was just an example)
Is it something similar?

[kerim371]

@mloubout you were right!

It seems there is typo in limit_model_to_receiver_area
I guess the snippet:

nx_min = Int(min_x ÷ model.d[1]) + 1

should be replaced by :

nx_min = Int((min_x-model.o[1]) ÷ model.d[1]) + 1

Do I need to create a PR?

I tried with this modification and it seems the modeling goes fine and consumes way less memory:

[mloubout]

You can open a PR yes but make sure you are on the latest JUDI version for that.