diff --git a/makefile.in b/makefile.in index 213a38d4..7f6332f5 100644 --- a/makefile.in +++ b/makefile.in @@ -38,7 +38,7 @@ speedupflag=yes # use compiler optimizations # # -blasflag=yes +blasflag=no sacflag=no @@ -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))) diff --git a/src/Acquisition.jl b/src/Acquisition.jl index 7ecb3d25..8d058ec2 100644 --- a/src/Acquisition.jl +++ b/src/Acquisition.jl @@ -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. diff --git a/src/Fdtd.jl b/src/Fdtd.jl index 155caefa..b0e8a973 100644 --- a/src/Fdtd.jl +++ b/src/Fdtd.jl @@ -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 diff --git a/src/Interpolation.jl b/src/Interpolation.jl index a332cf44..919a73d1 100644 --- a/src/Interpolation.jl +++ b/src/Interpolation.jl @@ -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, @@ -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 @@ -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 @@ -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 diff --git a/src/Inversion.jl b/src/Inversion.jl index 0e04d719..c9b3b084 100644 --- a/src/Inversion.jl +++ b/src/Inversion.jl @@ -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); @@ -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 diff --git a/src/Models.jl b/src/Models.jl index 2e0801a9..e98a4605 100644 --- a/src/Models.jl +++ b/src/Models.jl @@ -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);