examples: switch to segy_scan for larger datasets

slimgroup · Jul 12, 2022 · d4d1aab · d4d1aab
1 parent 8ae4be8
commit d4d1aab
Show file tree

Hide file tree

Showing 2 changed files with 19 additions and 8 deletions.
diff --git a/examples/scripts/lsrtm_2D.jl b/examples/scripts/lsrtm_2D.jl
@@ -19,7 +19,7 @@ model0 = Model((n[1], n[2]), (d[1], d[2]), (o[1], o[2]), m0)
 if ~isfile("$(JUDI.JUDI_DATA)/marmousi_2D.segy")
     ftp_data("ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/2DLSRTM/marmousi_2D.segy")
 end
-block = segy_read("$(JUDI.JUDI_DATA)/marmousi_2D.segy")
+block = segy_scan(JUDI.JUDI_DATA, "marmousi_2D.segy", ["GroupX","GroupY","RecGroupElevation","SourceSurfaceElevation","dt"])
 d_lin = judiVector(block)   # linearized observed data
 
 # Set up wavelet
@@ -41,17 +41,25 @@ Mr = judiTopmute(model0.n, 52, 10)
 
 #' set up number of iterations
 niter = parse(Int, get(ENV, "NITER", "10"))
+# Default to 64, 5 for CI only with NITER=1
+nsrc = 5 * parse(Int, get(ENV, "NITER", "$(q.nsrc ÷ 5)"))
+indsrc = randperm(q.nsrc)[1:nsrc]
 lsqr_sol = zeros(Float32, prod(n))
 
 # LSQR
-Ml = judiMarineTopmute2D(30, d_lin.geometry)
-lsqr!(lsqr_sol, Ml*J*Mr, Ml*d_lin; maxiter=niter)
+dinv = d_lin[indsrc]
+Jinv = J[indsrc]
+
+Ml = judiMarineTopmute2D(30, dinv.geometry)
+lsqr!(lsqr_sol, Ml*Jinv*Mr, Ml*dinv; maxiter=niter)
 
 # Save final velocity model, function value and history
 h5open("lsrtm_marmousi_lsqr_result.h5", "w") do file
     write(file, "x", reshape(lsqr_sol, model0.n))
 end
 
 # Plot final image
-figure(); imshow(copy(adjoint(reshape(lsqr_sol, model0.n))), extent = (0, 7.99, 3.19, 0), cmap = "gray", vmin = -3e-2, vmax = 3e-2)
-title("LS-RTM with LSQR"); xlabel("Lateral position [km]"); ylabel("Depth [km]")
+figure(); imshow(reshape(lsqr_sol, model0.n)', extent = (0, 7.99, 3.19, 0), cmap = "gray", vmin = -3e-2, vmax = 3e-2)
+title("LS-RTM with LSQR")
+xlabel("Lateral position [km]")
+ylabel("Depth [km]")
diff --git a/examples/scripts/splsrtm_2D.jl b/examples/scripts/splsrtm_2D.jl
@@ -19,7 +19,7 @@ model0 = Model((n[1], n[2]), (d[1], d[2]), (o[1], o[2]), m0)
 if ~isfile("$(JUDI.JUDI_DATA)/marmousi_2D.segy")
     ftp_data("ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/2DLSRTM/marmousi_2D.segy")
 end
-block = segy_read("$(JUDI.JUDI_DATA)/marmousi_2D.segy")
+block = segy_scan(JUDI.JUDI_DATA, "marmousi_2D.segy", ["GroupX","GroupY","RecGroupElevation","SourceSurfaceElevation","dt"])
 d_lin = judiVector(block)   # linearized observed data
 
 # Set up wavelet
@@ -72,5 +72,8 @@ h5open("lsrtm_marmousi_breg_result.h5", "w") do file
 end
 
 # Plot final image
-figure(); imshow(copy(adjoint(reshape(solb.x, model0.n))), extent = (0, 7.99, 3.19, 0), cmap = "gray", vmin = -3e-2, vmax = 3e-2)
-title("SPLS-RTM with Linearized Bregman"); xlabel("Lateral position [km]"); ylabel("Depth [km]")
+figure()
+imshow(reshape(solb.x, model0.n)', extent = (0, 7.99, 3.19, 0), cmap = "gray", vmin = -3e-2, vmax = 3e-2)
+title("SPLS-RTM with Linearized Bregman")
+xlabel("Lateral position [km]")
+ylabel("Depth [km]")