fix surface gather

Project.toml

@@ -1,14 +1,14 @@
 name = "ImageGather"
 uuid = "355d8124-6b2e-49b5-aab5-cdbc0a5fccbe"
 authors = ["mloubout <mathias.louboutin@gmail.com>"]
-version = "0.2.6"
+version = "0.2.7"
 
 [deps]
 JUDI = "f3b833dc-6b2e-5b9c-b940-873ed6319979"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
-JUDI = "3.3.1"
+JUDI = "3.3.9"
 julia = "1"
 
 [extras]

examples/layers_cig.jl

@@ -15,8 +15,8 @@ for i=1:12
     v[:,25*i+1:end] .= 1.5f0 + i*.25f0
 end
 v0 = imfilter(v, Kernel.gaussian(5))
-v0_low = .95f0 .* imfilter(v, Kernel.gaussian(45))
-v0_high = 1.05f0 .* imfilter(v, Kernel.gaussian(45))
+v0_low = .95f0 .* imfilter(v, Kernel.gaussian(5))
+v0_high = 1.05f0 .* imfilter(v, Kernel.gaussian(5))
 
 # Slowness squared [s^2/km^2]
 m = (1f0 ./ v).^2

src/implementation.py

@@ -45,7 +45,7 @@ def cig_grad(model, src_coords, wavelet, rec_coords, res, offsets, ic="as", spac
     hs = tuple(h.shape[0] for h in ohs.values())
     hd = tuple(h.indices[0] for h in ohs.values())
     gradm = Function(name="gradm", grid=model.grid, shape=(*hs, *u.shape[1:]),
-                    dimensions=(*hd, *model.grid.dimensions))
+                     dimensions=(*hd, *model.grid.dimensions))
 
     uh, vh = shifted_wf(u, v, ohs)
     g_expr = grad_expr(gradm, uh, vh, model, ic=ic)

src/surface_gather.jl

@@ -1,4 +1,4 @@
-import JUDI: AbstractModel
+import JUDI: AbstractModel, rlock_pycall, devito
 
 export surface_gather, double_rtm_cig
 
@@ -42,6 +42,8 @@ Compute the single shot contribution to the surface offset gather via double rtm
 
 """
 function double_rtm_cig(model_full, q::judiVector, data::judiVector, offs, options)
+    GC.gc(true)
+    devito.clear_cache()
     # Load full geometry for out-of-core geometry containers
     data.geometry = Geometry(data.geometry)
     q.geometry = Geometry(q.geometry)
@@ -56,11 +58,11 @@ function double_rtm_cig(model_full, q::judiVector, data::judiVector, offs, optio
 
     # Set up Python model
     modelPy = devito_model(model, options)
-    dtComp = get_dt(model; dt=options.dt_comp)
+    dtComp = convert(Float32, modelPy."critical_dt")
 
     # Extrapolate input data to computational grid
-    qIn = time_resample(make_input(q), q.geometry, dtComp)[1]
-    res = time_resample(make_input(data), data.geometry, dtComp)[1]
+    qIn = time_resample(make_input(q), q.geometry, dtComp)
+    res = time_resample(make_input(data), data.geometry, dtComp)
 
     # Set up coordinates
     src_coords = setup_grid(q.geometry, model.n)  # shifts source coordinates by origin
@@ -75,10 +77,10 @@ function double_rtm_cig(model_full, q::judiVector, data::judiVector, offs, optio
     res_o = res .* off_r'
     # Double rtm
 
-    rtm, rtmo, illum = JUDI.pylock() do
-        pycall(impl."double_rtm", Tuple{PyArray, PyArray, PyArray},
-               modelPy, qIn, src_coords, res, res_o, rec_coords, space_order=options.space_order)
-    end
+    rtm, rtmo, illum = rlock_pycall(impl."double_rtm", Tuple{PyArray, PyArray, PyArray},
+                                    modelPy, qIn, src_coords, res, res_o, rec_coords,
+                                    space_order=options.space_order)
+
     rtm = remove_padding(rtm, modelPy.padsizes)
     rtmo = remove_padding(rtmo, modelPy.padsizes)
     illum = remove_padding(illum, modelPy.padsizes)
@@ -94,6 +96,6 @@ function double_rtm_cig(model_full, q::judiVector, data::judiVector, offs, optio
     for (i, h) in enumerate(offs)
         soffs[:, :, i] .+= rtm .* delta_h(h_map, h, 2*diff(offs)[1])
     end
-    GC.gc()
+
     return soffs
 end