multicomponent receivers; forward and inversion of vector acoustic eq

src/FWI/FWI.jl

@@ -43,6 +43,11 @@ struct ModFdtdBorn end
 
 struct ModFdtdHBorn end
 
+# fields
+struct P end
+struct Vx end
+struct Vz end
+
 
 """
 FWI Parameters
@@ -132,21 +137,27 @@ include("xwfwi.jl")
 Convert the data `TD` to `Src` after time reversal.
 """
 function update_adjsrc!(adjsrc, δdat::Data.TD, adjacqgeom)
+	(adjsrc.fields != δdat.fields) && error("dissimilar fields")
 	nt=δdat.tgrid.nx
 	for i in 1:adjacqgeom.nss
-		for j in 1:length(δdat.fields)
+		for (j,field) in enumerate(δdat.fields)
 			wav=adjsrc.wav[i,j] 
 			dat=δdat.d[i,j]
 			for ir in 1:δdat.acqgeom.nr[i]
 				for it in 1:nt
-					@inbounds wav[it,ir]=dat[nt-it+1,ir]
+					@inbounds wav[it,ir]=value_adjsrc(dat[nt-it+1,ir],eval(field)())
 				end
 			end
 		end
 	end
 	return nothing
 end
 
+value_adjsrc(s, ::P) = s
+value_adjsrc(s, ::Vx) = -1.0*s # see adj tests
+value_adjsrc(s, ::Vz) = -1.0*s
+
+
 function generate_adjsrc(fields, tgrid, adjacqgeom)
 	adjsrc=Acquisition.Src_zeros(adjacqgeom, fields, tgrid)
 	return adjsrc
@@ -199,11 +210,11 @@ function Param(
 	       modm::Models.Seismic;
 	       # other optional 
 	       tgrid_obs::Grid.M1D=deepcopy(tgrid),
-	       recv_fields=[:P],
+	       rfields=[:P],
 	       igrid::Grid.M2D=nothing,
 	       igrid_interp_scheme::Symbol=:B2,
 	       mprecon_factor::Float64=1.0,
-	       dobs::Data.TD=Data.TD_zeros(recv_fields,tgrid_obs,acqgeom),
+	       dobs::Data.TD=Data.TD_zeros(rfields,tgrid_obs,acqgeom),
 	       dprecon=nothing,
 	       tlagssf_fracs=0.0,
 	       tlagrf_fracs=0.0,
@@ -260,13 +271,14 @@ function Param(
 	adjacqgeom = AdjGeom(acqgeom)
 
 	# generate adjoint sources
-	adjsrc=generate_adjsrc(recv_fields, tgrid, adjacqgeom)
+	adjsrc=generate_adjsrc(rfields, tgrid, adjacqgeom)
 
 	# generating forward and adjoint modelling engines
 	# to generate modelled data, border values, etc.
 	# most of the parameters given to this are dummy
 	paf=Fdtd.Param(npw=2, model=modm, born_flag=true,
 		acqgeom=[acqgeom, adjacqgeom], acqsrc=[acqsrc, adjsrc], 
+		rfields=rfields,
 		sflags=[3, 2], rflags=[1, 1],
 		backprop_flag=1, 
 		tgridmod=tgrid, gmodel_flag=true, verbose=verbose, illum_flag=false, nworker=nworker)
@@ -286,17 +298,17 @@ function Param(
 	end
 
 	# create Parameters for data misfit
-	#coup=Coupling.TD_delta(dobs.tgrid, tlagssf_fracs, tlagrf_fracs, recv_fields, acqgeom)
+	#coup=Coupling.TD_delta(dobs.tgrid, tlagssf_fracs, tlagrf_fracs, rfields, acqgeom)
 	# choose the data misfit
 	#	(iszero(tlagssf_fracs)) 
  	# tlagssf_fracs==[0.0] | tlagssf_fracs=[0.0])
-	# paTD=Data.P_misfit(Data.TD_zeros(recv_fields,tgrid,acqgeom),dobs,w=dprecon,coup=coup, func_attrib=optims[1]);
+	# paTD=Data.P_misfit(Data.TD_zeros(rfields,tgrid,acqgeom),dobs,w=dprecon,coup=coup, func_attrib=optims[1]);
 
 	if(iszero(dobs))
 		Random.randn!(dobs) # dummy dobs, update later
 	end
 
-	paTD=Data.P_misfit(Data.TD_zeros(recv_fields,tgrid,acqgeom),dobs,w=dprecon);
+	paTD=Data.P_misfit(Data.TD_zeros(rfields,tgrid,acqgeom),dobs,w=dprecon);
 
 	paminterp=Interpolation.Kernel([modi.mgrid.x, modi.mgrid.z], [modm.mgrid.x, modm.mgrid.z], igrid_interp_scheme)
 
@@ -311,7 +323,7 @@ function Param(
 	     mx, mxm,
 	     paf,
 	     deepcopy(acqsrc), 
-	     Acquisition.Src_zeros(adjacqgeom, recv_fields, tgrid),
+	     Acquisition.Src_zeros(adjacqgeom, rfields, tgrid),
 	     deepcopy(acqgeom), 
 	     adjacqgeom,
 	     deepcopy(modm), deepcopy(modm0), modi, mod_initial,

src/FWI/core.jl

@@ -186,7 +186,7 @@ function Fadj!(pa::Param)
 	pa.paf.c.illum_flag=false
 
 	# force boundaries in first pw and back propagation for second pw
-	pa.paf.c.sflags=[3,2] 
+	pa.paf.c.sflags=[-2,2] 
 	pa.paf.c.backprop_flag=-1
 
 	# update source wavelets in paf using adjoint sources

src/Fdtd/Fdtd.jl

@@ -19,6 +19,12 @@ global const to = TimerOutput(); # create a timer object
 global const npml = 50
 global const nlayer_rand = 0
 
+# struct related to fields
+struct P end
+struct Vx end
+struct Vz end
+
+# 
 #As forward modeling method, the 
 #finite-difference method is employed. 
 #It uses a discrete version of the two-dimensional isotropic acoustic wave equation.
@@ -64,6 +70,7 @@ mutable struct Paramc
 	acqgeom::Vector{Acquisition.Geom}
 	acqsrc::Vector{Acquisition.Src}
 	abs_trbl::Vector{Symbol}
+	sfields::Vector{Vector{Symbol}}
 	isfields::Vector{Vector{Int64}}
 	sflags::Vector{Int64} 
 	irfields::Vector{Int64}
@@ -361,13 +368,21 @@ function Param(;
 
 	# necessary that nss and fields should be same for all nprop
 	nss = acqgeom[1].nss;
-	sfields = [acqsrc[ipw].fields for ipw in 1:npw]
-	isfields = [findall(in([:P, :Vx, :Vz]), acqsrc[ipw].fields) for ipw in 1:npw]
+	sfields=[acqsrc[ipw].fields for ipw in 1:npw]
+	isfields = [Array{Int}(undef,length(sfields[ipw])) for ipw in 1:npw]
+	for ipw in 1:npw
+		for (i,field) in enumerate(sfields[ipw])
+			isfields[ipw][i]=findfirst(isequal(field), [:P,:Vx,:Vz])
+		end
+	end
 	fill(nss, npw) != [getfield(acqgeom[ip],:nss) for ip=1:npw] ? error("different supersources") : nothing
 
 	# create acquisition geometry with each source shooting 
 	# at every unique receiver position
-	irfields = findall(in([:P, :Vx, :Vz]), rfields)
+	irfields = Array{Int}(undef,length(rfields))
+	for (i,field) in enumerate(rfields)
+		irfields[i]=findfirst(isequal(field), [:P,:Vx,:Vz])
+	end
 
 
 	#(verbose) &&	println(string("\t> number of super sources:\t",nss))	
@@ -457,7 +472,7 @@ function Param(;
 	activepw=[ipw for ipw in 1:npw]
 	pac=Paramc(jobname,npw,activepw,
 	    exmodel,model,
-	    acqgeom,acqsrc,abs_trbl,isfields,sflags,
+	    acqgeom,acqsrc,abs_trbl,sfields,isfields,sflags,
 	    irfields,rflags,δt,δxI,δzI,
             nx,nz,nt,δtI,δx24I,δz24I,a_x,b_x,k_xI,a_x_half,b_x_half,k_x_halfI,a_z,b_z,k_zI,a_z_half,b_z_half,k_z_halfI,
 	    modtt, modttI,modrr,modrrvx,modrrvz,

src/Fdtd/source.jl

@@ -20,7 +20,7 @@ function fill_wavelets!(iss::Int64, wavelets::Array{Array{Float64,2},2}, acqsrc:
 					wavelets[ipw,it][is,ifield] = source_term
 				end
 			elseif(sflags[ipw] == -1)
-				"use this to add source sink fro previous case "
+				"source sink for 1"
 				"ϕ[t] = -s[t]"
 				for it=2:snt
 					source_term = acqsrc[ipw].wav[iss,ifield][it,is]
@@ -45,12 +45,10 @@ function fill_wavelets!(iss::Int64, wavelets::Array{Array{Float64,2},2}, acqsrc:
 				if(nt > snt)
 					wavelets[ipw,snt+1:end][is,ifield] = wavelets[ipw,snt][is,ifield]
 				end
-			elseif(sflags[ipw] == 3)
-				"use this to add source sink: need during adjoint propagation from boundary"
-				"multiplication with -1 for subtraction"
-				"time reversal"
-				"as the source wavelet has to be subtracted before the propagation step, I shift here by one sample"
-				"ϕ[t] = "
+			elseif(sflags[ipw] == -2)
+				"source sink for 2"
+				# multiplication with -1 for subtraction #time reversal
+				# as the source wavelet has to be subtracted before the propagation step, I shift here by one sample"
 				source_term_stack = 0.0;
 				for it=1:snt-1
 					source_term_stack += (acqsrc[ipw].wav[iss,ifield][it,is] .* δt)
@@ -80,48 +78,42 @@ struct Source_B0 end
 	isz2=pass[issp].isz2
 	ssprayw=pass[issp].ssprayw
 	modttI=pac.modttI
+	modrr=pac.modrrvz
 	"""
 	adding source to pressure field at [it] 
 	"""
 	for ipw in pac.activepw
+		sfields=pac.sfields[ipw]
+		isfields=pac.isfields[ipw]
 	if(pac.sflags[ipw] ≠ 0) # only if sflags is non-zero
 		pw=p[ipw]	
-		for (ifields, ifield) in enumerate(pac.isfields[ipw])
+		for (iff, ifield) in enumerate(isfields)
 		@simd for is = 1:acqgeom[ipw].ns[iss]
 			"""
 			use wavelets at [it], i.e., sum of source terms
 			until [it-1]
 			division of source term with δx and δz (see Jan's fdelmodc manual)
 			"""
-			source_term = wavelets[ipw,it][is, ifields] * pac.δt * pac.δxI * pac.δzI
+			source_term = wavelets[ipw,it][is, iff] * pac.δt * pac.δxI * pac.δzI
 
 			"""
 			multiplication with modttI
 			"""
 			pw[isz1[ipw][is], isx1[ipw][is],ifield] += 
-				source_term * 
-				ssprayw[ipw][1,is] * 
-				modttI[isz1[ipw][is], isx1[ipw][is]]  
+				source(source_term,ssprayw[ipw][1,is], pac, isz1[ipw][is], isx1[ipw][is],eval(sfields[iff])())
 			pw[isz1[ipw][is], isx2[ipw][is],ifield] += 
-				source_term * 
-				ssprayw[ipw][2,is] * 
-				modttI[isz1[ipw][is], isx2[ipw][is]]
+				source(source_term,ssprayw[ipw][2,is], pac, isz1[ipw][is], isx2[ipw][is],eval(sfields[iff])())
 			pw[isz2[ipw][is], isx1[ipw][is],ifield] += 
-				source_term * 
-				ssprayw[ipw][3,is] * 
-				modttI[isz2[ipw][is], isx1[ipw][is]]
+				source(source_term,ssprayw[ipw][3,is], pac, isz2[ipw][is], isx1[ipw][is],eval(sfields[iff])())
 			pw[isz2[ipw][is], isx2[ipw][is],ifield] += 
-				source_term * 
-				ssprayw[ipw][4,is] * 
-				modttI[isz2[ipw][is], isx2[ipw][is]]
+				source(source_term,ssprayw[ipw][4,is], pac, isz2[ipw][is], isx2[ipw][is],eval(sfields[iff])())
 		end
 		end
 	end
 	end
 end
 
 
-
 # This routine ABSOLUTELY should not allocate any memory, called inside time loop.
 @inbounds @fastmath function add_source!(it::Int64, issp::Int64, iss::Int64, pac::Paramc, pass::Vector{Paramss}, pap::Paramp, ::Source_B0)
 	# aliases
@@ -136,25 +128,38 @@ end
 	adding source to pressure field at [it] 
 	"""
 	for ipw in pac.activepw
+		sfields=pac.sfields[ipw]
+		isfields=pac.isfields[ipw]
 	if(pac.sflags[ipw] ≠ 0) # only if sflags is non-zero
 		pw=p[ipw]
-		for (ifields, ifield) in enumerate(pac.isfields[ipw])
+		for (iff, ifield) in enumerate(isfields)
 		@simd for is = 1:acqgeom[ipw].ns[iss]
 			"""
 			use wavelets at [it], i.e., sum of source terms
 			until [it-1]
 			division of source term with δx and δz (see Jan's fdelmodc manual)
 			"""
-			source_term = wavelets[ipw,it][is, ifields] * pac.δt * pac.δxI * pac.δzI
+			source_term = wavelets[ipw,it][is, iff] * pac.δt * pac.δxI * pac.δzI
 
 			"""
 			multiplication with modttI
 			"""
-			pw[isz1[ipw][is], isx1[ipw][is],ifield] += source_term * modttI[isz1[ipw][is], isx1[ipw][is]]  
+			pw[isz1[ipw][is], isx1[ipw][is],ifield] += source(source_term,1.0,pac,isz1[ipw][is],isx1[ipw][is],eval(sfields[iff])())
 		end
 	end
 	end
 	end
 end
 
 
+# on pressure grid
+source(source_term, spray, pac, iz, ix, ::P) = source_term * spray * pac.modttI[iz,ix]
+
+# on Vx grid
+source(source_term, spray, pac, iz, ix, ::Vx) = source_term * spray * pac.modrrvx[iz,ix]
+
+# on Vz grid
+source(source_term, spray, pac, iz, ix, ::Vz) = source_term * spray * pac.modrrvz[iz,ix]
+
+
+

src/Gallery/fwi.jl

@@ -1,7 +1,7 @@
 
 # create pizza problem
 
-function xfwi_problem(attrib::Symbol; born_flag=false)
+function xfwi_problem(attrib::Symbol; born_flag=false, rfields=[:Vx,:Vz,:P])
 
 	if(attrib==:pizza)
 		# starting model
@@ -33,6 +33,7 @@ function xfwi_problem(attrib::Symbol; born_flag=false)
 
 	if(born_flag)
 		pa = FWI.Param(acqsrc, acqgeom, tgrid, FWI.ModFdtdBorn(),
+		       rfields=rfields,
 		       model0, modm0=model0, modm_obs=model, 
 		       igrid_interp_scheme=igrid_interp_scheme, 
 		       igrid=igrid, parameterization=parameterization, verbose=false);
@@ -41,6 +42,7 @@ function xfwi_problem(attrib::Symbol; born_flag=false)
 	else
 		pa = FWI.Param(acqsrc, acqgeom, tgrid, FWI.ModFdtd(),
 		       model0, modm_obs=model, modm0=model0, 
+		       rfields=rfields,
 		       igrid_interp_scheme=igrid_interp_scheme, 
 		       igrid=igrid, parameterization=parameterization, verbose=false);
 	end

test/FWI/born_map.jl

@@ -1,43 +1,45 @@
 
-pa, model=JG.xfwi_problem(:pizza, born_flag=true)
+for rfields in [[:P], [:Vx], [:Vz]]
+	pa, model=JG.xfwi_problem(:pizza, born_flag=true, rfields=rfields)
 
 
-F=JF.operator_Born(pa);
+	F=JF.operator_Born(pa);
 
-x1=randn(size(F,2)) 
-x2=randn(size(F,2)) 
-x12=x1.+x2
+	x1=randn(size(F,2)) 
+	x2=randn(size(F,2)) 
+	x12=x1.+x2
 
 
-d12=F*x12
-d1=F*x1
-δmodtt1=copy(pa.paf.c.δmodtt)
-d2=F*x2
-δmodtt2=copy(pa.paf.c.δmodtt)
+	d12=F*x12
+	d1=F*x1
+	δmodtt1=copy(pa.paf.c.δmodtt)
+	d2=F*x2
+	δmodtt2=copy(pa.paf.c.δmodtt)
 
 
-d12new=d1.+d2
+	d12new=d1.+d2
 
-f=Misfits.error_squared_euclidean!(nothing, d12, d12new, nothing, norm_flag=true)
+	f=Misfits.error_squared_euclidean!(nothing, d12, d12new, nothing, norm_flag=true)
 
-@test f<1e-30
+	@test f<1e-29
 
 
-function adjtest()
-	x=randn(size(F,2))
-	y=randn(size(F,1))
-	a=LinearAlgebra.dot(y,F*x)
-	b=LinearAlgebra.dot(x,adjoint(F)*y)
-	c=LinearAlgebra.dot(x, transpose(F)*F*x)
-	println("adjoint test: ", a, "\t", b)       
-	@test isapprox(a,b,rtol=1e-6)
-	println("must be positive: ", c)
-	@test c>0.0
-end
+	function adjtest()
+		x=randn(size(F,2))
+		y=randn(size(F,1))
+		a=LinearAlgebra.dot(y,F*x)
+		b=LinearAlgebra.dot(x,adjoint(F)*y)
+		c=LinearAlgebra.dot(x, transpose(F)*F*x)
+		println("adjoint test: ", a, "\t", b)       
+		@test isapprox(a,b,rtol=1e-6)
+		println("must be positive: ", c)
+		@test c>0.0
+	end
 
 
-for i in 1:5
-	adjtest()
+	for i in 1:3
+		adjtest()
+	end
 end