Skip to content

Commit

Permalink
faster interpolation routines
Browse files Browse the repository at this point in the history
  • Loading branch information
pawbz committed Jul 13, 2017
1 parent 035cbc2 commit 197ba5c
Show file tree
Hide file tree
Showing 6 changed files with 82 additions and 43 deletions.
5 changes: 3 additions & 2 deletions makefile.in
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ speedupflag=yes
# use compiler optimizations
#
#
blasflag=yes
blasflag=no

sacflag=no

Expand All @@ -47,7 +47,8 @@ staticflag=no
doubleflag=yes

#fftlib = -L/scratch/intelics/mkl/lib/intel64 -lfftw3_threads -lfftw3
fftlib = -lfftw3 -lfftw3_threads
#fftlib = -lfftw3 -lfftw3_threads
fftlib =
#fftlib = -L/math/home/pawbz/fftw-3.3.5/lib -lfftw3_threads -lfftw3
#fftlib = -L/math/home/pawbz/fftw-3.3.5/lib -lfftw3
ifeq ($(strip $(sacflag)),$(strip $(yes)))
Expand Down
24 changes: 24 additions & 0 deletions src/Acquisition.jl
Original file line number Diff line number Diff line change
Expand Up @@ -181,6 +181,30 @@ function Geom_boundary(acqgeom::Geom,
end


"""
Check if all the sources and receivers in `Geom` are within the model
# Return
* `true` if all the positions within the model, `false` otherwise
"""
function Geom_check(geom::Geom, mgrid::Grid.M2D)
xmin, zmin, xmax, zmax = mgrid.x[1], mgrid.z[1], mgrid.x[end], mgrid.z[end]


checkvec=fill(false, geom.nss)
for iss=1:geom.nss
checkvec[iss] = any(
vcat(
(((geom.sx[iss]-xmin).*(xmax-geom.sx[iss])) .< 0.0),
(((geom.sz[iss]-zmin).*(zmax-geom.sz[iss])) .< 0.0),
(((geom.rx[iss]-xmin).*(xmax-geom.rx[iss])) .< 0.0),
(((geom.rz[iss]-zmin).*(zmax-geom.rz[iss])) .< 0.0),
) )
end
return !(any(checkvec))
end

"""
A fixed spread acquisition has same set of sources and
receivers for each supersource.
Expand Down
5 changes: 4 additions & 1 deletion src/Fdtd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -152,11 +152,14 @@ Author: Pawan Bharadwaj
# return
#endif

# all the propagating wavefield should have same sources, receivers, ? check that?
# all the propagating wavefields should have same supersources? check that?

# check dimension of model and model_pert

# check if all sources are receivers are inside model
any(![Acquisition.Geom_check(acqgeom[ip], model.mgrid) for ip=1:npropwav]) ? error("sources or receivers not inside model") : nothing



length(acqgeom) != npropwav ? error("acqgeom size") : nothing
length(acqsrc) != npropwav ? error("acqsrc size") : nothing
Expand Down
79 changes: 42 additions & 37 deletions src/Interpolation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,40 +12,45 @@ Return n indices in order
Cannot find a julia method which does, this.
If a faster method is found, replace it later.
"""
#function indminn(x::AbstractVector{Float64}, val::Float64, n::Int64=1)
# # using enumerate to avoid indexing
function indminn(x::AbstractVector{Float64}, val::Float64, n::Int64=1)
# using enumerate to avoid indexing
ivec = fill(0,n);
xc = zeros(n);
for inn in 1:n
min_i = 0
min_x = Inf
for (i, xi) in enumerate(x)
dist = abs(xi - val)
if ((dist < min_x))
min_x = dist
min_i = i
xc[]
end
end
xc[inn] = x[min_i]
x[min_i] = typemax(Float64)
ivec[inn] = min_i
end
x[ivec] = xc
return sort(ivec)
end

## slower version for large vectors
#function indminn(x::AbstractVector{Float64}, val::Float64, n::Int64)
# xa =abs(x-val)
# ((n >= 1) & (n <= length(x))) ? nothing : error("invalid n")
# ivec = [];
# for inn in 1:n
# min_i = 0
# min_x = Inf
# for (i, xi) in enumerate(x)
# dist = abs(xi - val)
# if ((dist < min_x) & (!(i in ivec)))
# min_x = dist
# min_i = i
# end
# end
# push!(ivec, min_i)
# xc = [];
# for i in 1:n
# ii = indmin(xa);
# push!(ivec, ii)
# push!(xc, x[ii])
# xa[ii] = typemax(Float64);
# end
# xa[ivec] = xc
# return sort(ivec)
#end

# slower version for large vectors
function indminn(x::AbstractVector{Float64}, val::Float64, n::Int64)
xa =abs(x-val)
((n >= 1) & (n <= length(x))) ? nothing : error("invalid n")
ivec = [];
xc = [];
for i in 1:n
ii = indmin(xa);
push!(ivec, ii)
push!(xc, x[ii])
xa[ii] = typemax(Float64);
end
xa[ivec] = xc
return sort(ivec)
end

function interp_spray!(xin::Vector{Float64}, yin::Vector{Float64},
xout::Vector{Float64}, yout::Vector{Float64},
attrib::Symbol,
Expand All @@ -66,12 +71,12 @@ function interp_spray!(xin::Vector{Float64}, yin::Vector{Float64},

yout[:] = zero(Float64);
if(attrib == :interp)
for i in eachindex(xout)
@simd for i in eachindex(xout)
ivec=indminn(xin,xout[i], np);
interp_func(view(xin, ivec), view(yin, ivec), view(xout, i), view(yout,i))
end
elseif(attrib == :spray)
for i in eachindex(xin)
@simd for i in eachindex(xin)
ivec=indminn(xout,xin[i], np);
spray_func(view(xout, ivec), view(yout,ivec),view(xin, i), view(yin, i))
end
Expand Down Expand Up @@ -99,15 +104,15 @@ function interp_spray!(xin::Array{Float64,1}, zin::Array{Float64,1}, yin::Array{
if(attrib == :interp)
y_x=zeros(Float64, length(zin), length(xout))
# first along x
for iz in eachindex(zin)
for ix in eachindex(xout)
for ix in eachindex(xout)
@simd for iz in eachindex(zin)
ivec=indminn(xin,xout[ix], np);
interp_func(view(xin, ivec), view(yin, iz,ivec), view(xout,ix), view(y_x,iz,ix))
end
end
# then along z
for ix in eachindex(xout)
for iz in eachindex(zout)
@simd for iz in eachindex(zout)
ivec=indminn(zin,zout[iz], np);
interp_func(view(zin, ivec), view(y_x, ivec,ix), view(zout,iz), view(yout,iz,ix))
end
Expand All @@ -116,14 +121,14 @@ function interp_spray!(xin::Array{Float64,1}, zin::Array{Float64,1}, yin::Array{
y_z = zeros(Float64, length(zout), length(xin))
# first along z
for ix in eachindex(xin)
for iz in eachindex(zin)
@simd for iz in eachindex(zin)
ivec=indminn(zout,zin[iz], np);
spray_func(view(zout, ivec), view(y_z,ivec,ix), view(zin ,iz), view(yin, iz,ix))
end
end
# then along x
for iz in eachindex(zout)
for ix in eachindex(xin)
for ix in eachindex(xin)
@simd for iz in eachindex(zout)
ivec=indminn(xout,xin[ix], np);
spray_func(view(xout, ivec), view(yout,iz,ivec), view(xin, ix), view(y_z, iz,ix))
end
Expand Down
6 changes: 6 additions & 0 deletions src/Inversion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -680,6 +680,9 @@ function Seismic_gx!(gmodm::Models.Seismic,
end
end

"""
Convert the data `TD` to `Src` after time reversal.
"""
function AdjSrc(δdat::Data.TD,
)
adjacq = AdjGeom(δdat.acqgeom);
Expand All @@ -699,6 +702,9 @@ function AdjGeom(geomin::Acquisition.Geom)
geomout.sx = geomin.rx; geomout.sz = geomin.rz;
geomout.ns = geomin.nr;

geomout.rx = geomin.sx; geomout.rz = geomin.sz;
geomout.nr = geomin.ns;

return geomout
end

Expand Down
6 changes: 3 additions & 3 deletions src/Models.jl
Original file line number Diff line number Diff line change
Expand Up @@ -418,11 +418,11 @@ function Seismic_interp_spray!(mod::Seismic, mod_out::Seismic, attrib::Symbol, B

"loop over fields in `Seismic`"
Interpolation.interp_spray!(mod.mgrid.x, mod.mgrid.z, mod.χvp,
mod_out.mgrid.x, mod_out.mgrid.z, mod_out.χvp, attrib)
mod_out.mgrid.x, mod_out.mgrid.z, mod_out.χvp, attrib, Battrib)
Interpolation.interp_spray!(mod.mgrid.x, mod.mgrid.z, mod.χvs,
mod_out.mgrid.x, mod_out.mgrid.z, mod_out.χvs, attrib)
mod_out.mgrid.x, mod_out.mgrid.z, mod_out.χvs, attrib, Battrib)
Interpolation.interp_spray!(mod.mgrid.x, mod.mgrid.z, mod.χρ,
mod_out.mgrid.x, mod_out.mgrid.z, mod_out.χρ, attrib)
mod_out.mgrid.x, mod_out.mgrid.z, mod_out.χρ, attrib, Battrib)

mod_out.vp0 = copy(mod.vp0); mod_out.vs0 = copy(mod.vs0);
mod_out.ρ0 = copy(mod.ρ0);
Expand Down

0 comments on commit 197ba5c

Please sign in to comment.