diff --git a/.github/workflows/runtests.yml b/.github/workflows/runtests.yml index 095fef1..24a431f 100644 --- a/.github/workflows/runtests.yml +++ b/.github/workflows/runtests.yml @@ -22,7 +22,7 @@ jobs: fail-fast: false matrix: - version: ['1.6', '1.7', '1.8', '1.9', '1'] + version: ['1.6', '1.7', '1.8', '1.9', '1.10', '1'] os: [ubuntu-latest] steps: 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..8e6cdf4 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..edc22ca 100644 Binary files a/docs/img/cig_line.png and b/docs/img/cig_line.png differ diff --git a/examples/layers_cig.jl b/examples/layers_cig.jl index 94fc233..5208129 100644 --- a/examples/layers_cig.jl +++ b/examples/layers_cig.jl @@ -61,7 +61,7 @@ wavelet = ricker_wavelet(timeS, dtS, f0) q = judiVector(srcGeometry, wavelet) ################################################################################################### -opt = Options(space_order=16) +opt = Options(space_order=8) # Setup operators F = judiModeling(model, srcGeometry, recGeometry; options=opt) # Nonlinear modeling 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..1b5a42b 100644 --- a/src/implementation.py +++ b/src/implementation.py @@ -10,12 +10,17 @@ 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"): + +def double_rtm(model, wavelet, src_coords, res, res_o, rec_coords, ic="as", space_order=8): """ """ _, u, illum, _ = forward(model, src_coords, None, wavelet, illum=True, - save=True, t_sub=4) + save=True, t_sub=4, space_order=space_order) # RTMs rtm = gradient(model, res, rec_coords, u, ic=ic)[0] @@ -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..d8307dd 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,9 +71,8 @@ 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 - off_r = log.(abs.(data.geometry.xloc[1] .- q.geometry.xloc[1]) .+ scale) - inv_off(x) = exp.(x) .- scale + scale = 5f3 + off_r = abs.(data.geometry.xloc[1] .- q.geometry.xloc[1]) .+ scale # mute if mute @@ -82,16 +81,15 @@ 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, space_order=options.space_order) rtm = remove_padding(rtm, modelPy.padsizes) rtmo = remove_padding(rtmo, modelPy.padsizes) illum = remove_padding(illum, modelPy.padsizes) # offset map - h_map = inv_off(offset_map(rtm, rtmo)) + h_map = offset_map(rtm, rtmo; scale=scale) rtm = laplacian(rtm) rtm[illum .> 0] ./= illum[illum .> 0]