How to generate plane wave using JUDI

[tjk9501]

Hello!

I want to use the JUDI software to setup a source which generates plane waves, but it seems that the default sources defined in JUDI are all point sources, can anyone tells me how to setup a plane wave source in JUDI?

Thanks!

[mloubout]

If you want to generate a plane wave, the easiest way is to use judiWeights with ones on the plane you want the plane wave at. See e.g for usage of judi weights:
https://github.com/slimgroup/JUDI.jl/blob/master/examples/scripts/modeling_extended_source_2D.jl

[tjk9501]

Could you help me see why the above code cannot output the wavefield data on the entire computation grid? It seems that using codes in the examples can only output wavefields at receiver positions.

[mloubout]

I'll have a look it should be working and tested

[tjk9501]

Hello!
Have you figured out using which expression can output wavefields on the entire 120 by 100 2D simulation domain?
Here are all the Packages we used on the ubuntu platform:
julia: v1.6.6
JUDI: v2.6.8
devito: v4.6.2
python: v3.7.13

[tjk9501]

If we query the type of d_sim, it will output the following message:

julia>d_sim
(JOLI.joLinearOperator{Float32,Float32}, (6024000, 1), "(forward wave equation*judiExtendedSource)")

[mloubout]

Could you try a newer version if JUDI? 2.6 is fairly old

[tjk9501]

Hello!
I tried newer version of JUDI, and here are the software setups on my computer:

Julia: v1.8.5
JUDI: v3.2.1

But when I retrive the wavefields data using

julia> Pw = judiLRWF(dt, wavelet)

julia prompts the following error:

Error showing value of type judiLRWF{Float32}:
ERROR: type judiWavelet has no field name
Stacktrace:
  [1] getproperty(P::judiLRWF{Float32}, s::Symbol)
    @ JUDI ~/.julia/packages/JUDI/wJSeg/src/TimeModeling/LinearOperators/lazy.jl:121
  [2] display(A::judiLRWF{Float32})
    @ JOLI ~/.julia/packages/JOLI/4APU1/src/joAbstractLinearOperator/base_functions.jl:24
  [3] #invokelatest#2
    @ ./essentials.jl:729 [inlined]
  [4] invokelatest
    @ ./essentials.jl:726 [inlined]
  [5] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
    @ REPL ~/Program/julia-1.8.5/share/julia/stdlib/v1.8/REPL/src/REPL.jl:296
  [6] (::REPL.var"#45#46"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL ~/Program/julia-1.8.5/share/julia/stdlib/v1.8/REPL/src/REPL.jl:278
  [7] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL ~/Program/julia-1.8.5/share/julia/stdlib/v1.8/REPL/src/REPL.jl:521
  [8] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL ~/Program/julia-1.8.5/share/julia/stdlib/v1.8/REPL/src/REPL.jl:276
  [9] (::REPL.var"#do_respond#66"{Bool, Bool, REPL.var"#77#87"{REPL.LineEditREPL, REPL.REPLHistoryProvider}, REPL.LineEditREPL, REPL.LineEdit.Prompt})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL ~/Program/julia-1.8.5/share/julia/stdlib/v1.8/REPL/src/REPL.jl:857
 [10] #invokelatest#2
    @ ./essentials.jl:729 [inlined]
 [11] invokelatest
    @ ./essentials.jl:726 [inlined]
 [12] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit ~/Program/julia-1.8.5/share/julia/stdlib/v1.8/REPL/src/LineEdit.jl:2510
 [13] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL ~/Program/julia-1.8.5/share/julia/stdlib/v1.8/REPL/src/REPL.jl:1248
 [14] (::REPL.var"#49#54"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL ./task.jl:484

I use the code example ''extended_source_lsqr.jl"
and the code snippet is

# Example for basic 2D modeling:
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, [email protected]
# Date: January 2017
#

using JUDI, SegyIO, LinearAlgebra, PyPlot, IterativeSolvers, JOLI

# 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.4f0
v[:,Int(round(end/2)):end] .= 5f0

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2

# Setup model structure
nsrc = 1	# number of sources
model = Model(n, d, o, m)

# 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
time = 2000f0   # receiver recording time [ms]
dt = 4f0    # receiver sampling interval [ms]

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

# Source wavelet
f0 = 0.01f0     # kHz
wavelet = ricker_wavelet(time, dt, f0)

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

# Write shots as segy files to disk
opt = Options()

# Setup operators
Pr = judiProjection(recGeometry)
F = judiModeling(model; options=opt)

# Random weights (size of the model)
w = judiWeights(randn(Float32, model.n))

# Create operator for injecting the weights, multiplied by the provided wavelet(s)
Pw = judiLRWF(dt, wavelet)

# Model observed data w/ extended source
lambda = 1f2
I = joDirac(prod(model.n), DDT=Float32, RDT=Float32)
F = Pr*F*adjoint(Pw)
F̄ = [F; lambda*I]

# Simultaneous observed data
d_sim = F*w

[tjk9501]

I tried your instructions, and julia stop prompt the Pw operator error, however, if I run the following code

F = judiModeling(model; options=opt)
Pw = judiLRWF(dt, wavelet);
F = F*adjoint(Pw)
d_sim = F*w

JUDI prompts the following error(at 'd_sim = F*w' expression)

julia> d_sim = F*w
ERROR: MethodError: no method matching devito_interface(::PyCall.PyObject, ::Matrix{Float32}, ::Matrix{Float32}, ::Nothing, ::Nothing, ::Nothing, ::JUDIOptions, ::Bool, ::Bool)
Closest candidates are:
  devito_interface(::PyCall.PyObject, ::Geometry, ::Array, ::Nothing, ::Nothing, ::Nothing, ::JUDIOptions, ::Bool, ::Bool) at ~/.julia/packages/JUDI/wJSeg/src/TimeModeling/Modeling/python_interface.jl:76
  devito_interface(::PyCall.PyObject, ::Nothing, ::Array, ::Nothing, ::Nothing, ::Nothing, ::JUDIOptions, ::Bool, ::Bool) at ~/.julia/packages/JUDI/wJSeg/src/TimeModeling/Modeling/python_interface.jl:103
  devito_interface(::PyCall.PyObject, ::Array, ::Array, ::Geometry, ::Nothing, ::Nothing, ::JUDIOptions, ::Bool, ::Bool) at ~/.julia/packages/JUDI/wJSeg/src/TimeModeling/Modeling/python_interface.jl:153
  ...
Stacktrace:
 [1] time_modeling(model_full::Model, srcGeometry::Matrix{Float32}, srcData::Matrix{Float32}, recGeometry::Nothing, recData::Nothing, dm::Nothing, op::Symbol, options::JUDIOptions, fw::Bool)
   @ JUDI ~/.julia/packages/JUDI/wJSeg/src/TimeModeling/Modeling/time_modeling_serial.jl:39
 [2] propagate(F::judiPointSourceModeling{Float32, :forward}, q::judiWeights{Float32})
   @ JUDI ~/.julia/packages/JUDI/wJSeg/src/TimeModeling/Modeling/propagation.jl:9
 [3] run_and_reduce(func::Function, #unused#::Nothing, nsrc::Int64, arg_func::JUDI.var"#234#235"{judiPointSourceModeling{Float32, :forward}, judiWeights{Float32}})
   @ JUDI ~/.julia/packages/JUDI/wJSeg/src/TimeModeling/Modeling/propagation.jl:32
 [4] multi_src_propagate(F::judiPointSourceModeling{Float32, :forward}, q::judiWeights{Float32})
   @ JUDI ~/.julia/packages/JUDI/wJSeg/src/TimeModeling/Modeling/propagation.jl:67
 [5] *(F::judiPointSourceModeling{Float32, :forward}, q::judiWeights{Float32})
   @ JUDI ~/.julia/packages/JUDI/wJSeg/src/TimeModeling/LinearOperators/operators.jl:173
 [6] top-level scope
   @ REPL[33]:1

[mloubout]

Ok thanks, looking into it, appreciate the patience

[mloubout]

Hum, looks like this case (output wavefield from extended source) isn't implemented. I will be working on adding it, sorry for the above I thought it was already there.

[mloubout]

Would you mind opening an issue so I can keep track of it?

[tjk9501]

Ok

[tjk9501]

Hello!
Before the update of newer versions, I tried the following method to record wavefields: In the following code I deploy receivers on each grid point of the simulation domain, but executing the code results in process being killed by julia. I don't know whether it is because the simulation exceeds memory upper bounds set in Julia (because similar setting can be run in other finite difference programs), and I want to find a way to make the code running. Here is the code snippet

using JUDI, SegyIO, LinearAlgebra, Plots
# 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.4f0
# v[:,Int(round(end/2)):end] .= 3f0
# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
nsrc = 2        # number of sources
model = Model(n, d, o, m)
nxrec = n[1]
nsrc = 1

#Setup Receiver Positions On Each Grid Point (Record Wavefields on the entire 2D simulation domain)
xrec = range(0f0, stop=1190f0, length=nxrec)
yrec = 0f0
zrec = range(0f0, stop=0f0, length=nxrec)

for index=1:n[2]-1
	global xrec = [xrec;xrec]
	ztemp = range(index*d[2],stop=index*d[2],length=n[1])
	global zrec = [zrec;ztemp]
end
#End of receiver setups

time = 1000f0   # receiver recording time [ms]
dt = 4f0    # receiver sampling interval [ms]
# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dt, t=time, nsrc=nsrc)
# Source wavelet
f0 = 0.01f0     # kHz
wavelet = ricker_wavelet(time, dt, f0)
xsrc = range(60f0,stop=60f0,length=nsrc)
ysrc = 0f0
zsrc = range(60f0,stop=60f0,length=nsrc)
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dt, t=time)
q = judiVector(srcGeometry,wavelet)
opt = Options(save_data_to_disk = true, return_array=false, dt_comp=1.0, free_surface=true, file_path="/home/tjk/Desktop",file_name = "SeismicTest")
nt = 251
Pr = judiProjection(recGeometry)
F = judiModeling(model; options=opt)
Ps = judiProjection(srcGeometry)
F = Pr*F*adjoint(Ps)
d_sim = F*q

[mloubout]

In the case of a point source (judiProjection * judiVecvtor) luke in this case getting a wavefield out is supported so you can do u = F * adjoint(Ps) * q. And if you want a plane wave you can set the source geometry as all the points along a plane (i.e fix ysrc/zrc and set xsrc = 0f0:d[1]:(n[1]-1)*d[1])

The options save_data_to_disk is however only for shot records (saved to SegY) so you can't use it for wavefields. See for example
https://github.com/slimgroup/JUDI.jl/blob/master/examples/scripts/modeling_wavefields_2D.jl

[tjk9501]

Thanks! Here is the case, I want to use the extended source function, but I want to verify how many actuating points are sufficient to generate an approximate plane wave. So, in practice, I want to use extended source and deploy receivers everywhere to check the wavefields. In this example I use single source to test whether this receiver deploying method is applicable.

[tjk9501]

Thanks! I have found there is a bug in defining the receivers, fixing the bug will make the code run smoothly and get the desired wavefields.