diff --git a/Project.toml b/Project.toml index 034769e..4cfbf6d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,14 +1,14 @@ name = "ImageGather" uuid = "355d8124-6b2e-49b5-aab5-cdbc0a5fccbe" authors = ["mloubout "] -version = "0.3.0" +version = "0.4.0" [deps] JUDI = "f3b833dc-6b2e-5b9c-b940-873ed6319979" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [compat] -JUDI = "3.4.6" +JUDI = "4" julia = "1" [extras] diff --git a/docs/img/cig_cdp.png b/docs/img/cig_cdp.png index 3a46cdd..023ad9f 100644 Binary files a/docs/img/cig_cdp.png and b/docs/img/cig_cdp.png differ diff --git a/docs/img/cig_line.png b/docs/img/cig_line.png index 7f92828..d4776fa 100644 Binary files a/docs/img/cig_line.png and b/docs/img/cig_line.png differ diff --git a/src/ImageGather.jl b/src/ImageGather.jl index 4460cfa..ec0332d 100644 --- a/src/ImageGather.jl +++ b/src/ImageGather.jl @@ -1,17 +1,20 @@ module ImageGather using JUDI - using JUDI.DSP, JUDI.PyCall + using JUDI.DSP, JUDI.PythonCall import Base: getindex, * import JUDI: judiAbstractJacobian, judiMultiSourceVector, judiComposedPropagator, judiJacobian, make_input, propagate import JUDI.LinearAlgebra: adjoint - const impl = PyNULL() + const impl = PythonCall.pynew() + + IGPath = pathof(ImageGather) function __init__() - pushfirst!(PyVector(pyimport("sys")."path"),dirname(pathof(ImageGather))) - copy!(impl, pyimport("implementation")) + pyimport("sys").path.append(dirname(IGPath)) + PythonCall.pycopy!(impl, pyimport("implementation")) + set_devito_config("autopadding", false) end # Utility functions include("utils.jl") diff --git a/src/implementation.py b/src/implementation.py index f4614ec..f370933 100644 --- a/src/implementation.py +++ b/src/implementation.py @@ -10,6 +10,11 @@ from kernels import wave_kernel from utils import opt_op +try: + from devitopro import * # noqa +except ImportError: + pass + def double_rtm(model, wavelet, src_coords, res, res_o, rec_coords, ic="as"): """ @@ -24,12 +29,12 @@ def double_rtm(model, wavelet, src_coords, res, res_o, rec_coords, ic="as"): def cig_grad(model, src_coords, wavelet, rec_coords, res, offsets, ic="as", - space_order=8, dims=None, illum=False): + space_order=8, dims=None, illum=False, t_sub=1): """ """ so = max(space_order, np.max(np.abs(offsets)) // model.grid.spacing[0]) - u = forward(model, src_coords, None, wavelet, - space_order=(space_order, so, so), save=True)[1] + _, u, _, _ = forward(model, src_coords, None, wavelet, t_sub=t_sub, + illum=illum, space_order=(space_order, so, so), save=True) # Setting adjoint wavefield v = wavefield(model, space_order, fw=False, tfull=True) diff --git a/src/subsurface_gather.jl b/src/subsurface_gather.jl index 1ea9a31..cdc416a 100644 --- a/src/subsurface_gather.jl +++ b/src/subsurface_gather.jl @@ -82,7 +82,7 @@ function propagate(J::judiExtendedJacobian{T, :born, O}, q::AbstractArray{T}, il modelPy = devito_model(J.model, J.options) nh = [length(J.offsets) for _=1:length(J.dims)] dmd = reshape(dm, nh..., J.model.n...) - dtComp = convert(Float32, modelPy."critical_dt") + dtComp = pyconvert(Float32, modelPy.critical_dt) # Extrapolate input data to computational grid qIn = time_resample(srcData, srcGeometry, dtComp) @@ -92,9 +92,9 @@ function propagate(J::judiExtendedJacobian{T, :born, O}, q::AbstractArray{T}, il rec_coords = setup_grid(recGeometry, J.model.n) # shifts rec coordinates by origin # Devito interface - dD = JUDI.wrapcall_data(impl."cig_lin", modelPy, src_coords, qIn, rec_coords, - dmd, J.offsets, ic=J.options.IC, space_order=J.options.space_order, dims=J.dims) - dD = time_resample(dD, dtComp, recGeometry) + dD = impl.cig_lin(modelPy, src_coords, qIn, rec_coords, dmd, J.offsets, + ic=J.options.IC, space_order=J.options.space_order, dims=J.dims) + dD = time_resample(PyArray(dD), dtComp, recGeometry) # Output shot record as judiVector return judiVector{Float32, Matrix{Float32}}(1, recGeometry, [dD]) end @@ -107,7 +107,7 @@ function propagate(J::judiExtendedJacobian{T, :adjoint_born, O}, q::AbstractArra # Set up Python model modelPy = devito_model(J.model, J.options) - dtComp = convert(Float32, modelPy."critical_dt") + dtComp = pyconvert(Float32, modelPy.critical_dt) # Extrapolate input data to computational grid qIn = time_resample(srcData, srcGeometry, dtComp) @@ -118,11 +118,10 @@ function propagate(J::judiExtendedJacobian{T, :adjoint_born, O}, q::AbstractArra rec_coords = setup_grid(recGeometry, J.model.n) # shifts rec coordinates by origin # Devito - g = JUDI.pylock() do - pycall(impl."cig_grad", PyArray, modelPy, src_coords, qIn, rec_coords, dObserved, J.offsets, - illum=false, ic=J.options.IC, space_order=J.options.space_order, dims=J.dims) - end - g = remove_padding_cig(g, modelPy.padsizes; true_adjoint=J.options.sum_padding) + g = impl.cig_grad(modelPy, src_coords, qIn, rec_coords, dObserved, J.offsets, + illum=false, ic=J.options.IC, space_order=J.options.space_order, dims=J.dims, + t_sub=J.options.subsampling_factor) + g = remove_padding_cig(PyArray(g), pyconvert(Tuple, modelPy.padsizes); true_adjoint=J.options.sum_padding) return g end diff --git a/src/surface_gather.jl b/src/surface_gather.jl index d895569..db5badf 100644 --- a/src/surface_gather.jl +++ b/src/surface_gather.jl @@ -1,4 +1,4 @@ -import JUDI: AbstractModel, rlock_pycall, devito +import JUDI: AbstractModel, devito, wrapcall_data export surface_gather, double_rtm_cig @@ -60,7 +60,7 @@ function double_rtm_cig(model_full, q::judiVector, data::judiVector, offs, optio # Set up Python model modelPy = devito_model(model, options) - dtComp = convert(Float32, modelPy."critical_dt") + dtComp = pyconvert(Float32, modelPy.critical_dt) # Extrapolate input data to computational grid qIn = time_resample(make_input(q), q.geometry, dtComp) @@ -71,7 +71,7 @@ function double_rtm_cig(model_full, q::judiVector, data::judiVector, offs, optio rec_coords = setup_grid(data.geometry, model.n) # shifts rec coordinates by origin # Src-rec offsets - scale = 1f1 + scale = 1f2 off_r = log.(abs.(data.geometry.xloc[1] .- q.geometry.xloc[1]) .+ scale) inv_off(x) = exp.(x) .- scale @@ -82,9 +82,8 @@ function double_rtm_cig(model_full, q::judiVector, data::judiVector, offs, optio res_o = res .* off_r' # Double rtm - rtm, rtmo, illum = rlock_pycall(impl."double_rtm", Tuple{PyArray, PyArray, PyArray}, - modelPy, qIn, src_coords, res, res_o, rec_coords, - ic=options.IC) + rtm, rtmo, illum = wrapcall_data(impl."double_rtm", modelPy, qIn, src_coords, res, res_o, rec_coords, + ic=options.IC) rtm = remove_padding(rtm, modelPy.padsizes) rtmo = remove_padding(rtmo, modelPy.padsizes)