Need help with using SPG.

[zhangxiaoshuotttt]

Hello @mloubout @kerim371 ,

Please allow me to provide some details about the problems.

I am facing a strange issue while using SPG. Before using SPG, I attempted to use PGN and L-BFGS on my field data. As these methods worked fine, it seems that the field data is not the issue here.

However, when I tried to use SPG on field data, I encountered some errors. Before this, I tried the SPG on overthrust data, which worked well.

I have attached my scripts and the errors for your reference.

#####  Read data  #####
container = segy_scan(prestk_dir, prestk_file, ["SourceX", "SourceY", "GroupX", "GroupY", "RecGroupElevation", "SourceSurfaceElevation", "dt"])
d_obs = judiVector(container;segy_depth_key = "RecGroupElevation")

#####  Set Model  #####
model0 = Model(n,d,o,m0,rho=rho0, nb=300,abc_type="pml");

#####  Source function  ###### 
segy_depth_key_src = "SourceSurfaceElevation"

# Set up wavelet and source vector
src_geometry = Geometry(container; key = "source", segy_depth_key = segy_depth_key_src)

f0 = 0.005f0

rickerwavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], f0)

q = judiVector(src_geometry,rickerwavelet)
plot(q.data[1])

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

#####  Filter  ###### 
frq = 0.005 # kHz
# Bandpass filter
Ml_freq = judiFilter(d_obs.geometry, 0.001, frq*1000f0)
Mr_freq = judiFilter(q.geometry, 0.001, frq*1000f0)

############# JUDI options #############
########################################

global jopt = JUDI.Options(
    IC = "fwi",
    subsampling_factor=4, 
    #limit_m = true,
    #buffer_size = buffer_size,
    optimal_checkpointing=false,
    free_surface=true,  # free_surface is ON to model multiples as well
    space_order=16)     # increase space order for > 12 Hz source wavelet

# optimization parameters
niterations = 20
count = 0
fhistory = Vector{Float32}(undef, 0)
mute_reflections = false
mute_turning = false

#####  SETUP CONSTARAINTS ###### 

options=PARSDMM_options()
options.FL=Float32
options=default_PARSDMM_options(options,options.FL)
constraint = Vector{SetIntersectionProjection.set_definitions}()

# Bound constraints
vmin = 1.4
vmax = 5.2
vBoundCoef = 0.5
kr = 50

# Slowness squared [s^2/km^2]
mmin = (1f0 ./ vmax).^2
mmax = (1f0 ./ vmin).^2
mminArr = ones(Float32, size(model0)) .* mmin
mmaxArr = ones(Float32, size(model0)) .* mmax

############# Loss function #############
########################################
@everywhere function H(x)
    n = size(x, 1)
    σ = ifftshift(sign.(-n/2+1:n/2))
    y = imag(ifft(σ.*fft(x, 1), 1))
    return y
end

@everywhere envelope(x, y) = sum(abs2.((x - y) .+ 1im .* H(x - y)))
@everywhere denvelope(x, y) = Zygote.gradient(xs->envelope(xs, y), x)[1]
@everywhere myloss(x, y) = (envelope(x, y), denvelope(x, y))

@everywhere myloss(randn(Float32, 10, 10), randn(Float32, 10, 10))

############# Object function #############
########################################
function objective_function(m_update)
    global x, z, count, jopt, seabed_ind;
    count += 1
    print("Now is the round is: ",count,"\n")
    # Update model
    model0.m .= Float32.(reshape(m_update, size(model0)))
    if modeling_type == "bulk"
        model0.rho .= Float32.(reshape(rho_from_slowness(model0.m), size(model0)))
        model0.rho[:,1:seabed_ind] .= rhowater
    end

    # fwi function value and gradient
    indsrc = randperm(d_obs.nsrc)[1:batchsize]
    if mute_reflections
        fval, gradient = fwi_objective(model0, Mr_freq[indsrc]*q[indsrc], Ml_tur[indsrc]*Ml_freq[indsrc]*get_data(d_obs[indsrc]), options=jopt, misfit=myloss)
    elseif mute_turning
        fval, gradient = fwi_objective(model0, Mr_freq[indsrc]*q[indsrc], Ml_ref[indsrc]*Ml_freq[indsrc]*get_data(d_obs[indsrc]), options=jopt, misfit=myloss)
    else
        fval, grad = fwi_objective(model0, Mr_freq[indsrc]*q[indsrc], Ml_freq[indsrc]*get_data(d_obs[indsrc]), options=jopt, misfit=myloss)
    end
    gradient = grad.data
    #gradient = reshape(gradient, size(model0))
    gradient[:, 1:seabed_ind] .= 0f0
    gradient = .125f0*gradient/maximum(abs.(gradient))  # scale for line search

    push!(fhistory, fval)

    println("iteration: ", count, "\tfval: ", fval, "\tnorm: ", norm(gradient))

    return fval, gradient
end

############# Perform FWI #############
########################################
# Optimization parameters
batchsize = 1

count = 0

# Bound projection
proj(x) = reshape(median([vec(mminArr) vec(x) vec(mmaxArr)]; dims=2), size(model0))
Minput = vec(model0.m.data)
# FWI with SPG
options = spg_options(verbose=3, maxIter=niterations, memory=3)
sol = spg(objective_function, Minput, proj, options)

Here is the error

After going through the error, I checked the original code and suspected that a particular line may be causing the issue.

I compared the model data output to the field data output and attempted to match their types, but the error persisted.

[kerim371]

I will take a look a little bit later

[mloubout]

Can you try Minput = model0.m.data[:] . From the error log g::Vector, x:;Matrix the vec might not work properly

[kerim371]

Yes, I just tested and @mloubout found the the source of the problem. But @zhangxiaoshuotttt you only need to set:

Minput = model0.m.data

That way it works

[zhangxiaoshuotttt]

Hi @mloubout @kerim371 , I tried both of your advice but still have the same issue. By the way, I can get the gradient, but it seems like I can't get into the SPG process.

[mloubout]

What version of JUDI and SlimOptim are you on?

[zhangxiaoshuotttt]

Sorry, I forget about an important thing. The version is 3.3.8. By the way, to implement the PML condition, I use this version provided by Kerim. I use this version because the field data has a thin water layer(Only 80m). If I use version 3.3.10. Then, the simulated data has some waves reflected from boundaries.

[mloubout]

Few comments

• the version you are referring too is a proof of concept branch that isn't finished and doesn't implement different absorbing boundary.
• if there is issues with the latest JUDI it would be appreciated to open an issue about it with how to reproduce it so it can be fixed rather than keeping it for you.

Since both @kerim371 and I can run SPG fine I would make sure that your versiobs are up to date so that we can help

[zhangxiaoshuotttt]

I understand. I will try it immediately. Thank you very much.

[mloubout]

As of shallow water, I've been running very shallow water (<30m) field data for a few month with 3.3.10 and with the opened field data PR and it works completely fine so not sure what issue you are facing.

[zhangxiaoshuotttt]

Hi again，I have recently updated to JUDI 3.3.10 and have been running simulations with the updated code.

#### Parameters for model ######

n = (480, 300)   # (x,y,z) or (x,z)
d = (5, 5)
o = (0., 0.)
seabed=80
seabedpoint=convert(Int64,seabed/d[1])
# Velocity [km/s]
v = ones(Float32,n) .+ 0.5f0
v0 = ones(Float32,n) .+ 0.5f0
v[:,seabedpoint:150] .= 2.2f0
v[:,150:end] .= 2.5f0
rho = (v0 .+ .5f0) ./ 2;

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

# Setup model structure
model = Model(n, d, o, m)
model0 = Model(n, d, o, m0)

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

#### Parameters for Geometry ######
nxrec = 120
xrec = range(0f0, stop=(n[1]-1)*d[1], length=nxrec)
yrec = 0f0 # WE have to set the y coordiante to zero (or any number) for 2D modeling
zrec = range(80f0, step=0, length=nxrec)

# receiver sampling and recording time
timeD = 4000f0   # receiver recording time [ms]
dtD = 0.5f0    # receiver sampling interval [ms]
nsrc = 12
# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtD, t=timeD, nsrc=nsrc)

xsrc = convertToCell(range(35f0, stop=(n[1]-5)*d[1], length=nsrc))
ysrc = convertToCell(range(0f0, stop=0f0, length=nsrc))
zsrc = convertToCell(range(20f0, step=0, length=nsrc))

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

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

# setup wavelet
f0 = 0.01f0     # kHz
wavelet = ricker_wavelet(timeD, dtD, f0)
q = judiVector(srcGeometry, wavelet)

# Setup options
opt = Options(space_order=12,free_surface=true,limit_m=true,buffer_size = 20)

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

dobs = Pr*F0*adjoint(Ps)*q

plot_sdata(dobs[6].data[1],(0.0005,10),cmap="RdBu",name="Water 80 with multiples",save=OutDir)

I have also been examining my simulated shot gather.

I have attempted to modify the buffer_size parameter from 0 to 3e3, but it appears to not be working well. Specifically, I have concerns about the effectiveness of the absorbing boundary. As a result, my shot gather is not as good as I would like it to be. If the absorbing boundary works well, I may get the shot gather like this.

[mloubout]

5m grid for a 10Hz wavelet is quite overkill the boundary layer is definitely gonna be short for such a setup, I would either augment my frequency, coarse the grid or add more boundary points to make sure there is enough wavelength in the 40 point boundary layer

[zhangxiaoshuotttt]

I see. And I have tried it. It truly works. And it works. I forgot to et the nb in the model. This was a careless mistake. Thank you for your patience and time.

[zhangxiaoshuotttt]

Hi @mloubout @kerim371 ， I have tried spg with new JUDI 3.3.10 with this new option, and I don't change other parameter, It still give me the same error. It's wired.

# JUDI options
global jopt = JUDI.Options(
    IC = "fwi",
    limit_m = true,
    buffer_size = buffer_size,
    optimal_checkpointing=false,
    free_surface=true,  
    space_order=16)     

[kerim371]

Sorry I've been away from the devices.

The error mentionned above is most likely connected with the fact that some of the arguments are Vector and other are Matrix. The first thing I would test is to check that all the variables (model, grad, Minput and other) have appropriate types (i.e. typeof) and the size. Use reshape when it is needed. Is your model0.m is 1D array or it is 2D? Minput 1D or 2D?

I've tested your code and it was giving the error you mentionned in the beginning but then I've made changes that I said and the error disappeared.

Your JUDI.Options seems to be fine.

My setup: devito: 4.8.2 JUDI: 3.9 SlimOptim: 0.2.3

[zhangxiaoshuotttt]

Thank you for all of your advice @kerim371 @mloubout . I was able to identify the error - it was due to the type of my rho, which was Float64. To resolve the issue, I changed it to Float32. Once this was done, the code worked perfectly. I appreciate your patience and the time you took to assist me. Thanks again!

[mloubout]

Great, glad it got resolved. Could you open an issue so I don't forget there is a type mismatch that needs to be fixed.