slimgroup / JUDI.jl

Julia Devito inversion.
https://slimgroup.github.io/JUDI.jl
MIT License
97 stars 31 forks source link

judiVector_to_seisblock error #85

Closed ziyiyin97 closed 2 years ago

ziyiyin97 commented 2 years ago

This MFE includes a verbatim copy of https://github.com/slimgroup/JUDI.jl/blob/master/examples/scripts/modeling_basic_2D.jl to generate data. It errors when turning judiVector to SeisBlock

# Example for basic 2D modeling:
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using JUDI, SegyIO, LinearAlgebra

# Set up model structure
n = (120, 100)   # (x,y,z) or (x,z)
d = (10., 10.)
o = (0., 0.)

# Velocity [km/s]
v = ones(Float32,n) .+ 0.5f0
v0 = ones(Float32,n) .+ 0.5f0
v[:,Int(round(end/2)):end] .= 3.5f0
rho = (v0 .+ .5f0) ./ 2

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
m0 = (1f0 ./ v0).^2
dm = vec(m - m0)

# Setup info and model structure
nsrc = 2    # number of sources
model = Model(n, d, o, m)
model0 = Model(n, d, o, m0)

# Set up receiver geometry
nxrec = 120
xrec = range(50f0, stop=1150f0, length=nxrec)
yrec = 0f0
zrec = range(50f0, stop=50f0, length=nxrec)

# receiver sampling and recording time
timeR = 1000f0   # receiver recording time [ms]
dtR = 2f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

# Set up source geometry (cell array with source locations for each shot)
xsrc = convertToCell(range(400f0, stop=800f0, length=nsrc))
ysrc = convertToCell(range(0f0, stop=0f0, length=nsrc))
zsrc = convertToCell(range(200f0, stop=200f0, length=nsrc))

# source sampling and number of time steps
timeS = 1000f0  # ms
dtS = 2f0   # ms

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)

# setup wavelet
f0 = 0.01f0     # kHz
wavelet = ricker_wavelet(timeS, dtS, f0)
q = judiVector(srcGeometry, wavelet)

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry, recGeometry, model)
info = Info(prod(n), nsrc, ntComp)

###################################################################################################

# Write shots as segy files to disk
opt = Options(optimal_checkpointing=false, isic=false, subsampling_factor=2, dt_comp=1.0)

# Setup operators
Pr = judiProjection(info, recGeometry)
F = judiModeling(info, model; options=opt)
F0 = judiModeling(info, model0; options=opt)
Ps = judiProjection(info, srcGeometry)
J = judiJacobian(Pr*F0*adjoint(Ps), q)

# Nonlinear modeling
dobs = Pr*F*adjoint(Ps)*q

block = judiVector_to_SeisBlock(dobs, q; source_depth_key="SourceDepth")

the error is

ERROR: LoadError: MethodError: no method matching round(::Vector{Float32})
Closest candidates are:
  round(::Union{Dates.Day, Dates.Week, Dates.TimePeriod}, ::Union{Dates.Day, Dates.Week, Dates.TimePeriod}, ::RoundingMode{:NearestTiesUp}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/Dates/src/rounding.jl:267
  round(::Union{Dates.Day, Dates.Week, Dates.TimePeriod, Dates.TimeType}, ::Dates.Period) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/Dates/src/rounding.jl:282
  round(::Union{Dates.Day, Dates.Week, Dates.TimePeriod, Dates.TimeType}, ::Dates.Period, ::RoundingMode{:Down}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/Dates/src/rounding.jl:273
  ...
Stacktrace:
 [1] judiVector_to_SeisBlock(d::judiVector{Float32, Matrix{Float32}}, q::judiVector{Float32, Matrix{Float32}}; source_depth_key::String, receiver_depth_key::String)
   @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/Types/judiVector.jl:435
 [2] top-level scope
   @ ~/.julia/dev/JUDI/MFE.jl:79
 [3] include(fname::String)
   @ Base.MainInclude ./client.jl:444
 [4] top-level scope
   @ REPL[1]:1
in expression starting at /Users/francisyin/.julia/dev/JUDI/MFE.jl:79

and the problem really is yrec=0f0. If we turn it to yrec=range(0f0, stop=0f0, length=nxrec) then it works.

This in the end leads to further concern: should 2D seismic data has a yloc as [[0],[0]] (for 2 sources), or [[0,0,....0],[0,0,...,0]]? It's either the first one or the second one, and turning from judiVector to SeisBlock then back to judiVector must be an "identity" operation which remains the geometry.

Plan to fix this by #69

kerim371 commented 2 years ago

Hi,

Just in case I think the same error occur when I set opt = Options(;save_data_to_disk=true) in the same modeling_basic_2D.jl example

JUDI v2.6.4