updated a few more tests to v0.7

Project.toml

@@ -33,6 +33,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
 Signals = "6303bc30-01ab-52eb-8b10-60e47555a8d1"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/Born/Born.jl

@@ -4,6 +4,8 @@ using Grid
 import JuMIT.Acquisition
 import JuMIT.Models
 import JuMIT.Data
+using SpecialFunctions
+using FFTW
 
 mutable struct Param
 end
@@ -80,7 +82,7 @@ function mod(;
 			for it in 1:nt
 				wpow2[it]=complex(acqsrc.wav[iss,ifield][it,is])
 			end
-			fft!(wpow2) # source wavelet in the frequency domain
+			FFTW.fft!(wpow2) # source wavelet in the frequency domain
 
 			x = sx[is] - rx[ir]
 			z = sz[is] - rz[ir]
@@ -124,7 +126,7 @@ function mod(;
 		end
 
 		# back to time domain
-		ifft!(dpow2all)
+		FFTW.ifft!(dpow2all)
 
 		# truncate
 		for it in 1:nt

src/Gallery/Gallery.jl

@@ -107,7 +107,7 @@ function Seismic(attrib::Symbol, δ::Float64=0.0; verbose=false)
 		vp0=Models.bounds(vp,bfrac); 
 		vs0=Models.bounds(vs,bfrac); 
 		ρ0=Models.bounds(ρ, bfrac);
-		mgrid = Grid.M2D(0., 17000., 0., 3500.,size(vp,2),size(vp,1),40)
+		mgrid = Grid.M2D(0., 17000., 0., 3500.,size(vp,2),size(vp,1))
 		model= Models.Seismic(vp0, vs0, ρ0, Models.χ(vp,vp0,1), Models.χ(vs,vs0,1), Models.χ(ρ,ρ0,1), mgrid)
 	elseif(attrib == :seismic_marmousi2_high_res)
 		vp, h= IO.readsu(joinpath(marmousi_folder,"vp_marmousi-ii.su"))
@@ -117,7 +117,7 @@ function Seismic(attrib::Symbol, δ::Float64=0.0; verbose=false)
 		vp0=Models.bounds(vp,bfrac); 
 		vs0=Models.bounds(vs,bfrac); 
 		ρ0=Models.bounds(ρ, bfrac);
-		mgrid = Grid.M2D(0., 17000., 0., 3500.,size(vp,2),size(vp,1),40)
+		mgrid = Grid.M2D(0., 17000., 0., 3500.,size(vp,2),size(vp,1))
 		model= Models.Seismic(vp0, vs0, ρ0, Models.χ(vp,vp0,1), Models.χ(vs,vs0,1), Models.χ(ρ,ρ0,1), mgrid)
 
 	elseif(attrib == :seismic_marmousi2_xwell)

test/Fdtd/accuracy.jl

@@ -1,47 +1,3 @@
-model = J.Gallery.Seismic(:acou_homo2);
-model0 = J.Gallery.Seismic(:acou_homo2);
-J.Models.Seismic_addon!(model,randn_perc=1, fields=[:χvp,:χρ])
-acqgeom=J.Acquisition.Geom_fixed(model,5,10)
-acqsrc=J.Acquisition.Src_fixed_mod(acqgeom.nss,1,[:P],mod=model, nλ=3, tmaxfrac=1.0)
-tgrid=acqsrc.tgrid
-
-attrib_mod=JF.ModFdtdBorn()
-
-#for born_flag in [true, false]
-pa=JF.Param(acqsrc, acqgeom, tgrid, attrib_mod, model0, 
-modm_obs=model,  
-modm0=model0,
-igrid_interp_scheme=:B2,    
-igrid=Grid.M2D_resamp(model.mgrid, 300.,300.,),     
-parameterization=[:χvp, :χρ, :null],   verbose=false,
-nworker=1)
-
-
-pa_parallel=JF.Param(acqsrc, acqgeom, 
-	tgrid, attrib_mod, model0, 
-	modm_obs=model,  
-	modm0=model0,
-	igrid_interp_scheme=:B2,    
-	igrid=Grid.M2D_resamp(model.mgrid, 300.,300.,),     
-	parameterization=[:χvp, :χρ, :null],   verbose=false,
-	nworker=nothing)
-
-
-paerr=JuMIT.Data.P_misfit(pa_parallel.paTD.y, pa.paTD.y)
-err=JuMIT.Data.func_grad!(paerr)
-
-# normalize error
-error = err[1]/paerr.ynorm
-
-
-	result=JF.xfwi!(pa, JF.Migr())
-	result_parallel=JF.xfwi!(pa_parallel, JF.Migr())
-
-
-f=Misfits.error_squared_euclidean!(nothing, result[2], result_parallel[2], nothing, norm_flag=true)
-
-
-ffff
 
 model = JuMIT.Gallery.Seismic(:acou_homo1);
 acqgeom = JuMIT.Gallery.Geom(model.mgrid,:xwell);
@@ -53,7 +9,7 @@ acqsrc = JuMIT.Acquisition.Src_fixed(acqgeom.nss,1,[:P],wav,tgrid);
 
 vp0=mean(JuMIT.Models.χ(model.χvp,model.ref.vp,-1))
 ρ0=mean(JuMIT.Models.χ(model.χρ,model.ref.ρ,-1))
-rec1 = JuMIT.Analytic.mod(vp0=vp0,
+rec1 = JuMIT.Born.mod(vp0=vp0,
             model_pert=model,
 		     ρ0=ρ0,
 			 acqgeom=acqgeom, acqsrc=acqsrc, tgridmod=tgrid, src_flag=2)

test/Fdtd/backprop.jl

@@ -1,41 +1,44 @@
 
 model=JuMIT.Gallery.Seismic(:acou_homo1)
-JuMIT.Models.Seismic_addon!(model, randn_perc=0.001, fields=[:χvp,:χρ])
+JuMIT.Models.Seismic_addon!(model, randn_perc=0.1, fields=[:χvp,:χρ])
 #model.mgrid.npml=5;
 acqgeom =JuMIT.Acquisition.Geom_circ(nss=4,nr=100,rad=[990.,990.]);
 acqgeom =JuMIT.Acquisition.Geom_circ(nss=4,nr=100,rad=[0.,200.]);
 acqsrc=JuMIT.Acquisition.Src_fixed_mod(acqgeom.nss,1,[:P],mod=model, nλ=3, tmaxfrac=0.4)
-@time pa=JuMIT.Fdtd.Param(born_flag=false,npw=1, tgridmod=acqsrc.tgrid,
-#	abs_trbl=[:null],
-        gmodel_flag=false,
-		sflags=[1],
-        snaps_flag=true,
-	verbose=true,
-        backprop_flag=1,
-        illum_flag=true,acqgeom=[acqgeom], acqsrc=[acqsrc],
-        model=model);
-
-JuMIT.Fdtd.mod!(pa);
-rec1=deepcopy(pa.c.data[1])
-
-# change source flag and update wavelets in pa
-pa.c.sflags=[-1];
-JuMIT.Fdtd.update_acqsrc!(pa,[acqsrc])
-pa.c.backprop_flag=-1 # do backpropagation
-
-JuMIT.Fdtd.mod!(pa);
-rec2=deepcopy(pa.c.data[1])
-
-# time reverse
-JuMIT.Data.TD_tr!(rec2);
-
-# compare results
-# least-squares misfit
-paerr=JuMIT.Data.P_misfit(rec1, rec2)
-err = JuMIT.Data.func_grad!(paerr)
-
-# normalized error
-error = err[1]./paerr.ynorm
-
-# desired accuracy?
-@test error<1e-20
+
+for sflags in [[1,-1],[2,3]]
+	pa=JuMIT.Fdtd.Param(born_flag=false,npw=1, tgridmod=acqsrc.tgrid,
+	#	abs_trbl=[:null],
+		gmodel_flag=false,
+		sflags=[sflags[1]],
+		snaps_flag=true,
+		verbose=true,
+		backprop_flag=1,
+		illum_flag=true,acqgeom=[acqgeom], acqsrc=[acqsrc],
+		model=model);
+
+	JuMIT.Fdtd.mod!(pa);
+	rec1=deepcopy(pa.c.data[1])
+
+	# change source flag and update wavelets in pa
+	pa.c.sflags=[sflags[2]];
+	JuMIT.Fdtd.update_acqsrc!(pa,[acqsrc])
+	pa.c.backprop_flag=-1 # do backpropagation
+
+	JuMIT.Fdtd.mod!(pa);
+	rec2=deepcopy(pa.c.data[1])
+
+	# time reverse
+	JuMIT.Data.TD_tr!(rec2);
+
+	# compare results
+	# least-squares misfit
+	paerr=JuMIT.Data.P_misfit(rec1, rec2)
+	err = JuMIT.Data.func_grad!(paerr)
+
+	# normalized error
+	error = err[1]./paerr.ynorm
+
+	# desired accuracy?
+	@test error<1e-20
+end

test/Models/Models.jl

@@ -7,19 +7,19 @@ mod=randn(1000); mod0=2.0
 @btime JuMIT.Models.χg!(mod,mod0,-1)
 
 # some model
-mgrid = Grid.M2D(0.0, 10., 0., 10., .05,.05,2);
+mgrid = Grid.M2D(0.0, 10., 0., 10., .05,.05);
 model = JuMIT.Models.Seismic_zeros(mgrid);
 vp0=[2100.,2200.];vs0=[-1., -1.]; ρ0=[2100., 2300.]
 JuMIT.Models.adjust_bounds!(model, vp0,vs0,ρ0);
 
 JuMIT.Models.Seismic_addon!(model, randn_perc=1e-3)
 nznx = model.mgrid.nz * model.mgrid.nx;
 
-# test copy!
+# test copyto!
 model_new=deepcopy(model)
 JuMIT.Models.Seismic_addon!(model_new, randn_perc=1e-3)
 
-@btime copy!(model_new, model)
+@btime copyto!(model_new, model)
 @test isequal(model_new, model) 
 
 # allocation of model0
@@ -56,13 +56,3 @@ end
 
 
 
-
-# testing pad_trun
-a=randn(200,300)
-np=50
-b=zeros(200+2*np,300+2*np)
-@time JuMIT.Models.pml_pad_trun!(b, a,1);
-aa=similar(a);
-@time JuMIT.Models.pml_pad_trun!(b, aa,-1);
-
-@test aa==a

test/Models/param_adj.jl

@@ -1,6 +1,6 @@
 
 # some model
-mgrid = Grid.M2D(0.0, 10., 0., 10., .05,.05,2);
+mgrid = Grid.M2D(0.0, 10., 0., 10., .05,.05);
 model = JuMIT.Models.Seismic_zeros(mgrid);
 vp0=[2100.,2200.];vs0=[-1., -1.]; ρ0=[2100., 2300.]
 JuMIT.Models.adjust_bounds!(model, vp0,vs0,ρ0);

test/runtests.jl

@@ -1,5 +1,5 @@
 using Distributed
-#addprocs(5)
+addprocs(5)
 
 using LinearAlgebra
 using JuMIT
@@ -18,15 +18,15 @@ function initialize(fp)
 end
 
 folder="Models"
-for t in ["param_adj"]
+for t in ["Models", "param_adj"]
 	fp = joinpath(folder, string(t, ".jl"))
 	initialize(fp)
 end
 
 
 
 folder="Fdtd"
-for t in ["rho_projection"]
+for t in ["accuracy", "backprop", "rho_projection"]
 	fp = joinpath(folder, string(t, ".jl"))
 	initialize(fp)
 end