Fix #88, now skip oob recs since generated code skips then and output…

…s zeros
slimgroup · Jan 23, 2024 · 830829e · 830829e
1 parent eaec430
commit 830829e
Show file tree

Hide file tree

Showing 4 changed files with 4 additions and 7 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "JUDI"
 uuid = "f3b833dc-6b2e-5b9c-b940-873ed6319979"
 authors = ["Philipp Witte, Mathias Louboutin"]
-version = "3.3.9"
+version = "3.3.10"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

diff --git a/src/TimeModeling/Modeling/time_modeling_serial.jl b/src/TimeModeling/Modeling/time_modeling_serial.jl
@@ -36,11 +36,6 @@ function time_modeling(model_full::AbstractModel, srcGeometry::GeomOrNot, srcDat
         modelPy = devito_model(model, options, dm)
     end
 
-    # Remove receivers outside the modeling domain (otherwise leads to segmentation faults)
-    @juditime "Remove OOB src/rec" begin
-        recGeometry, recData = remove_out_of_bounds_receivers(recGeometry, recData, model)
-    end
-
     # Devito interface
     @juditime "Propagation" begin
         argout = devito_interface(modelPy, srcGeometry, srcData, recGeometry, recData, dm, options, illum, fw)

diff --git a/src/TimeModeling/Utils/auxiliaryFunctions.jl b/src/TimeModeling/Utils/auxiliaryFunctions.jl
@@ -242,6 +242,7 @@ function remove_out_of_bounds_receivers(recGeometry::Geometry{T}, recData::Matri
     if typeof(recGeometry.xloc[1]) <: Array
         idx_xrec = findall(x -> xmax >= x >= xmin, recGeometry.xloc[1])
         recGeometry.xloc[1] = recGeometry.xloc[1][idx_xrec]
+        length(recGeometry.yloc[1]) > 1 && (recGeometry.yloc[1] = recGeometry.yloc[1][idx_xrec])
         recGeometry.zloc[1] = recGeometry.zloc[1][idx_xrec]
         recData = recData[:, idx_xrec]
     end
@@ -250,6 +251,7 @@ function remove_out_of_bounds_receivers(recGeometry::Geometry{T}, recData::Matri
     if length(size(model)) == 3 && typeof(recGeometry.yloc[1]) <: Array
         ymin, ymax = origin(model, 2), origin(model, 2) + (size(model, 2) - 1)*spacing(model, 2)
         idx_yrec = findall(x -> ymax >= x >= ymin, recGeometry.yloc[1])
+        recGeometry.xloc[1] = recGeometry.xloc[1][idx_yrec]
         recGeometry.yloc[1] = recGeometry.yloc[1][idx_yrec]
         recGeometry.zloc[1] = recGeometry.zloc[1][idx_yrec]
         recData = recData[:, idx_yrec]

diff --git a/src/pysource/fields.py b/src/pysource/fields.py
@@ -129,7 +129,7 @@ def wavefield_subsampled(model, u, nt, t_sub, space_order=8):
     if t_sub > 1:
         time_subsampled = ConditionalDimension(name='t_sub', parent=model.grid.time_dim,
                                                factor=t_sub)
-        nsave = nt // t_sub + 1
+        nsave = int(np.ceil((nt + t_sub)/t_sub))
     else:
         return None
     wf_s = []