PML boundary condition

[kerim371]

Current implementation of damping boundary condition requires greately extending the computational area.

I hope PML do that better.

To incorporate PML abc we should probably have a list of damp func. As a first step I'm trying to leave current implementation of damp_op function but make Model.damp member as list. With that I currently get error: ValueError('Defaultdamp(x, y)is incompatible with other args asx_size=65, whilex_size=200is expected. Perhaps you forgot to overridedamp(x, y)?')

The example of PML in Devito

Appreciate any help

[kerim371]

I'm still struggling with this error. What is the possible reason of such error? I simply wrapped Model.damp member to a list of a single function. The full error output:

ERROR: PyError ($(Expr(:escape, :(ccall(#= /home/kerim/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'ValueError'>
ValueError('Default `damp(x, y)` is incompatible with other args as `x_size=65`, while `x_size=200` is expected. Perhaps you forgot to override `damp(x, y)`?')
  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 80, in forward
    summary = op(**kw)
  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/operator/operator.py", line 738, in __call__
    return self.apply(**kwargs)
  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/operator/operator.py", line 804, in apply
    args = self.arguments(**kwargs)
  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/operator/operator.py", line 645, in arguments
    args = self._prepare_arguments(**kwargs)
  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/operator/operator.py", line 549, in _prepare_arguments
    raise ValueError("Default `%s` is incompatible with other args as "

Stacktrace:
  [1] pyerr_check
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/exception.jl:75 [inlined]
  [2] pyerr_check
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/exception.jl:79 [inlined]
  [3] _handle_error(msg::String)
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/exception.jl:96
  [4] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/exception.jl:110 [inlined]
  [5] #107
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:43 [inlined]
  [6] disable_sigint
    @ ./c.jl:458 [inlined]
  [7] __pycall!
    @ ~/Documents/Colada/r/julia-1.6/.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}}, nargs::Int64, kw::PyCall.PyObject)
    @ PyCall ~/Documents/Colada/r/julia-1.6/.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}}, 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/ilqDX/src/pyfncall.jl:11
 [10] #pycall#112
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:80 [inlined]
 [11] #4
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:82 [inlined]
 [12] (::JUDI.var"#1#2"{JUDI.var"#4#5"{Tuple{PyCall.PyArray, PyCall.PyObject}, 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}}}})()
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:74
 [13] lock(f::JUDI.var"#1#2"{JUDI.var"#4#5"{Tuple{PyCall.PyArray, PyCall.PyObject}, 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}}}}, l::ReentrantLock)
    @ Base ./lock.jl:187
 [14] pylock(f::JUDI.var"#4#5"{Tuple{PyCall.PyArray, PyCall.PyObject}, 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}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:71
 [15] rlock_pycall(::PyCall.PyObject, ::Type{Tuple{PyCall.PyArray, 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/JUDI.jl:81
 [16] 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
 [17] 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:65
 [18] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/time_modeling_serial.jl:46 [inlined]
 [19] macro expansion
    @ ./timing.jl:287 [inlined]
 [20] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:141 [inlined]
 [21] time_modeling(model_full::JUDI.IsoModel{Float32, 2}, srcGeometry::GeometryIC{Float32}, srcData::Matrix{Float32}, recGeometry::GeometryIC{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:45
 [22] 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
 [23] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:36 [inlined]
 [24] macro expansion
    @ ./timing.jl:287 [inlined]
 [25] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:141 [inlined]
 [26] run_and_reduce(func::Function, #unused#::Nothing, nsrc::Int64, arg_func::JUDI.var"#267#268"{judiDataSourceModeling{Float32, :forward}, judiVector{Float32, Matrix{Float32}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:35
 [27] 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:85
 [28] *
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/LinearOperators/operators.jl:173 [inlined]
 [29] afoldl
    @ ./operators.jl:533 [inlined]
 [30] *(a::judiProjection{Float32}, b::judiModeling{Float32, :forward, JUDI.IsoModel{Float32, 2}}, c::JUDI.jAdjoint{judiProjection{Float32}}, xs::judiVector{Float32, Matrix{Float32}})
    @ Base ./operators.jl:560
 [31] top-level scope
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/examples/scripts/modeling_basic_2D.jl:118

[mloubout]

You need to make sure that

https://github.com/slimgroup/JUDI.jl/blob/83752f147e175608d98b2ffe2d4778d58281c1f1/src/pysource/models.py#L260

Contains your pml field or it will use the ones from EmptyModel which is a tiny model only used to create the operator

[kerim371]

@mloubout I ve added JUDIOptions.abc_type field with possible values damp and pml. Then Model have attributes like pmlx0,pmlx1,pmly0 etc if the JUDIOptions.abc_type == pml.

I'm not sure whether I correctly implemented Model.pml_op() method.

I have some troubles with implementing kernel equations with pml. If you could guide me I would like to work further on this PR.

[kerim371]

I think the abc_type should go in Model, like NHL, rather than the options

Done, abc_type now is a member of DiscreteGrid struct as well as nb.

[mloubout]

I have some troubles with implementing kernel equations with pml.

What's the main blocker?

[kerim371]

@mloubout basically I think it is better if I push my efforts with mistakes and then we discuss it. I will do that soon (I need few days).

But ultimately I can't understand the following: for acoustic media PML abc are defined using these three equations: Does that mean that in kernels.py I have to define three stencils (like stencil = solve(wmr * u.dt2 + damp * udt - ulaplace - q, u_n)) and pack them in list of size 3 pde = [Eq(u_n, stencil, subdomain=model.grid.subdomains['nofsdomain'])]?

And further I don't know yet what are PML equations for scalar/Q-factor/VTI/TTI models.

[mloubout]

yes it would return a list of three equations instead of one

[kerim371]

@mloubout hi,

From the Devito PML example I have to implement phi1/phi2 functions:

In Devito example they are defined as follows:

phi1 = TimeFunction(name="phi1",grid=grid,time_order=2,space_order=2,staggered=(x,z),dtype=np.float64)
phi2 = TimeFunction(name="phi2",grid=grid,time_order=2,space_order=2,staggered=(x,z),dtype=np.float64)

but as I see we don't set any values to them even if in the Devito examples it is written about these vars: "auxiliary variables, which will be non zero only in the absorption region"

Am I missing something or I can simply define them as TimeFunction(...) and then takes their derivatives etc?

[mloubout]

but as I see we don't set any values

Yes this is expected, these are just extra state variables for the PML region, similar to the wavefield 'uso there is no value set it's an extra PDE which at each time step updates the values of Phi based onu` and the pml equation.

[kerim371]

@mloubout some new obstacles appeared...

I define phi as Model attributes:

self.phi1 = TimeFunction(name="phi1",grid=self.grid,time_order=2,space_order=2,staggered=self.grid.dimensions)
self.phi2 = TimeFunction(name="phi2",grid=self.grid,time_order=2,space_order=2,staggered=self.grid.dimensions)

I have a feeling that I have to add them to Operator so they were compiled. If so what should be the right hand side?

After that I will have to compute something like phi1[t,x,z-1]: while I can get x,z from Grid.Dimensions, how can I retrieve time dimension t?

[mloubout]

have a feeling that I have to add them to [Operator]so they were compiled.

Not sure what you mean, these are just objects, the operator just generates the code for the equations you define with it You should have a look at:

https://nbviewer.org/github/devitocodes/devito/tree/master/examples/seismic/abc_methods/

For some example of oimplementation

[zhangxiaoshuotttt]

I tried this today. It works well for implementing damping boundary conditions. Thank you all.

[kerim371]

@zhangxiaoshuotttt hi,

What I did here is:

1. added abc_type argument to Model structure (possible values: damp and pml)
2. tried to add pml grid to src/pysource/models.py (I'm sure whether I did that correctly)

And yet to implement in src/pysource/kernels.py stencil for each type of modeling (acoustic, elastic etc) based on equations from Devito example.

That was pretty difficult for me so I stopped this PR for the now.