search

slimgroup / JUDI.jl

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

JUDI DepthScaling is in-place now #40

Closed ziyiyin97 closed 3 years ago

ziyiyin97 commented 3 years ago

JUDI Depth Scaling operator is in-place right now, which should not be the case! A minimal failing example is here

using JUDI.TimeModeling, LinearAlgebra, JOLI, Test, Printf
# Set up model structure
n = (401, 401)   # (x,y,z) or (x,z)
d = (2.5f0, 2.5f0)
o = (0f0, 0f0)

# Velocity [km/s]
v0 = 2f0*ones(Float32,n)
v = deepcopy(v0);
v[201,201] = 3f0;

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

# Setup info and model structure    # number of sources
model0 = Model(n, d, o, m0; nb = 200)

dm1 = deepcopy(dm)

S = judiDepthScaling(model0)

dm2 = S*dm

@test isapprox(dm1,dm)

The problem is at https://github.com/slimgroup/JUDI.jl/blob/869564526393262efb3fe2286d57d4ff25ce1a36/src/TimeModeling/seismic_preconditioners.jl#L194, it should be 

m = copy(reshape(m,model.n))

I am going to fix this with a PR and add a test case to it.