From 7705e8b164ed70341b69c9a51510904216fba617 Mon Sep 17 00:00:00 2001 From: Pawan Bharadwaj Date: Sat, 26 Oct 2024 16:54:13 +0530 Subject: [PATCH] Update version to v0.19.46 in marmousi-acoustic.jl notebook Remove unused import of ColorSchemes in GeoPhyInv.jl Increase _fd_npml and _fd_npextend values in GeoPhyInv.jl Add maxiter parameter to IterativeSolvers.lsmr! in fwi.jl Update dependencies in Project.toml Add tauii_buffer to P_x_worker_x_pw struct in fdtd.jl Update add_velocity_source! and add_stress_source! functions in propagate.jl Remove unused seriescolor line in plots.jl Update compute_stressii! function in advance_elastic.jl Update gradient_finitediff_testing.jl notebook --- Manifest.toml | 2 +- Project.toml | 1 - docs/make.jl | 1 - notebooks/gradient_finitediff_testing.jl | 24 +- notebooks/marmousi-acoustic.jl | 2 +- notebooks/reusing-seisforwexpt.jl | 52 ++- notebooks/sample_fwi.jl | 502 +++++++++++++++++++++++ src/GeoPhyInv.jl | 5 +- src/fdtd/advance_elastic.jl | 10 +- src/fdtd/diff2D.jl | 27 +- src/fdtd/diff3D.jl | 39 +- src/fdtd/fdtd.jl | 3 +- src/fdtd/propagate.jl | 4 +- src/fdtd/source.jl | 88 +++- src/fdtd/types.jl | 4 +- src/fwi/fwi.jl | 2 +- src/plots.jl | 5 +- 17 files changed, 721 insertions(+), 50 deletions(-) create mode 100644 notebooks/sample_fwi.jl diff --git a/Manifest.toml b/Manifest.toml index 16bc996..2ab1b5d 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.4" manifest_format = "2.0" -project_hash = "b93906f24b48daeaa115e6a8eb2d9b29862f8e80" +project_hash = "da2b1e181419a1d0edd7834abebe8c15fb72559a" [[deps.ASL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] diff --git a/Project.toml b/Project.toml index 247fdaf..fcb25b4 100644 --- a/Project.toml +++ b/Project.toml @@ -8,7 +8,6 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" -ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" diff --git a/docs/make.jl b/docs/make.jl index cb29d97..2435a74 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -3,7 +3,6 @@ ENV["GKSwstype"] = "100" using Documenter, GeoPhyInv, SparseArrays, LinearAlgebra, Random, LinearMaps using Plots, Test, Luxor -using ColorSchemes # standard setting for generating doc pages # @init_parallel_stencil(2, false, Float32, 2) diff --git a/notebooks/gradient_finitediff_testing.jl b/notebooks/gradient_finitediff_testing.jl index e85c87d..f9f074a 100644 --- a/notebooks/gradient_finitediff_testing.jl +++ b/notebooks/gradient_finitediff_testing.jl @@ -55,7 +55,7 @@ end # @time update!(m1, pa_mod, [:invK]) # ╔═╡ f49439df-3db3-48ae-bf9d-1387e3c14c8f -N = 2 +N = 100 # ╔═╡ 6d47ae1d-2554-4d42-af6c-02ce952bf14e @bind mpara Select([:invK, :rho]) @@ -180,7 +180,7 @@ isnan.(pa_inv.dobs[1].d[1]) begin g = similar(m) GeoPhyInv.gradient!(g, m, pa_inv) - heatmap(Array(pa_mod.c.gradients[1])) + heatmap(Array(reshape(g, N, N)), c=:seismic) end # ╔═╡ aa0f7f1b-65ff-4a0b-9fb0-e6660911fdc8 @@ -198,6 +198,9 @@ heatmap(Array(reshape(pa_inv.P' * m, Nz, Nx))) # ╔═╡ 01552458-5de7-4627-b417-21745db2ade1 heatmap(Array(reshape(m, 176, 176))) +# ╔═╡ 8ec00fb7-fe18-4741-b987-813ec97d73fe +heatmap(Array(reshape(gm, N, N)), c=:seismic) + # ╔═╡ cd5bff59-6466-43d7-bd10-795d6f55d4ef heatmap(Array(reshape(g, 2, 2)), c=:seismic) @@ -213,6 +216,9 @@ gKI = reshape(g, N, N) # ╔═╡ 630faaff-f7f0-46d7-a36f-31a54e589c8e heatmap(Array(gKI), c=:seismic) +# ╔═╡ c6007b1a-3819-4543-bab4-f0dc0c12e075 +gm1 = Array(reshape(gm, N, N)) + # ╔═╡ 839cf146-639c-4d13-9470-d25882837169 # prod(step.(pa_true.c.medium.grid)) # pa_true.c.srcwav[1][1].grid |> step |> abs2 @@ -223,6 +229,9 @@ GeoPhyInv.Fields(pa_mod.c.attrib_mod, "v") # ╔═╡ b14a5a36-94c4-4a6c-8b26-730919e3e5ae filter(x -> x ∉ Fields(pa_mod.c.attrib_mod, "d"), Fields(pa_mod.c.attrib_mod, "v")) +# ╔═╡ 60a9adc0-466d-4e54-9838-f01b0a2c50c5 +Array(gKI) ./ gm1 + # ╔═╡ 417388ad-8709-4d4a-9d81-d4bbb0351a92 heatmap(Array(pa_mod.c.gradients[1]), clim=(-1e-11, 1e-11), c=:seismic) @@ -238,16 +247,7 @@ end Js(GeoPhyInv.Data.Array(xs)) # ╔═╡ c6626e0f-4bf1-4e29-a0c3-1e7f23805f61 -gm = grad(central_fdm(2, 1, factor=1e8), Js, xs)[1] - -# ╔═╡ 8ec00fb7-fe18-4741-b987-813ec97d73fe -heatmap(Array(reshape(gm, N, N)), c=:seismic) - -# ╔═╡ c6007b1a-3819-4543-bab4-f0dc0c12e075 -gm1 = Array(reshape(gm, N, N)) - -# ╔═╡ 60a9adc0-466d-4e54-9838-f01b0a2c50c5 -Array(gKI) ./ gm1 +# gm = grad(central_fdm(2, 1, factor=1e8), Js, xs)[1] # ╔═╡ 9371b452-4ad6-4d07-ae98-c86239b6c110 gs = get_xs(g, xs) diff --git a/notebooks/marmousi-acoustic.jl b/notebooks/marmousi-acoustic.jl index 94fd26e..baed58f 100644 --- a/notebooks/marmousi-acoustic.jl +++ b/notebooks/marmousi-acoustic.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.19.21 +# v0.19.46 #> [frontmatter] #> title = "Marmousi" diff --git a/notebooks/reusing-seisforwexpt.jl b/notebooks/reusing-seisforwexpt.jl index be17fe1..f760e4d 100644 --- a/notebooks/reusing-seisforwexpt.jl +++ b/notebooks/reusing-seisforwexpt.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.19.21 +# v0.19.46 #> [frontmatter] #> title = "SeisForwExpt" @@ -8,14 +8,26 @@ using Markdown using InteractiveUtils +# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error). +macro bind(def, element) + quote + local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end + local el = $(esc(element)) + global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) + el + end +end + # ╔═╡ 5fc7d3bf-418d-4487-a2d6-bc07c736374b begin import Pkg Pkg.add("PlutoLinks") + Pkg.add("Revise") Pkg.add("PlutoUI") Pkg.add("PlutoTest") Pkg.add("Plots") Pkg.add("CUDA") + using Revise using PlutoLinks: @revise using PlutoUI, PlutoTest, Plots, CUDA end @@ -52,8 +64,10 @@ various bundles of medium parameters. We will first load a predefined medium, a medium = ElasticMedium(Homogeneous(), 5) # ╔═╡ ded21d43-99d6-453b-a858-e32f3f82e605 -# surface seismic acquisition with 3 supersources and 100 receivers -ageom = AGeom(medium.grid, :surf, SSrcs(1), Recs(100)) +begin + # surface seismic acquisition with 3 supersources and 100 receivers + ageom = AGeom(medium.grid, :xwell, SSrcs(1), Recs(100)); +end # ╔═╡ a47b1e5d-a08c-4750-b1f3-f8a62225a19f # lets choose a time grid @@ -84,14 +98,19 @@ md""" Now we have all the stuff necessary to allocate an instance of `SeisForwExpt`. """ +# ╔═╡ 14713ee7-6753-46a7-98b3-ed0b8ef48ae4 +tsnaps=tgrid + # ╔═╡ c7d96dda-17f5-4dbb-98bf-4c29fd784574 pa = SeisForwExpt( FdtdElastic(), + snaps_field=:vx, + tsnaps=tsnaps, medium=medium, ageom=ageom, srcwav=srcwav, tgrid=tgrid, - rfields=[:vz], + rfields=[:vz, :vx], verbose=true, ); @@ -103,7 +122,7 @@ Before modeling, we will now generate a perturbed medium, which is "similar" to # ╔═╡ df4a4072-a91a-471a-88d9-b74ed58e7823 begin medium_box = deepcopy(medium) - update!(medium_box, [:vp, :rho, :vs], rectangle=[[-500, -500], [500, 500]], perc=5.0) # perturbed velocity box + update!(medium_box, [:vp, :rho, :vs], rectangle=[[-500, -500], [500, 500]], perc=20.0) # perturbed velocity box end @@ -127,6 +146,12 @@ begin data_box = deepcopy(pa[:data]) # copy data end +# ╔═╡ e9c4cb8b-20bf-43c3-8da5-9ecb2e0ddfd2 +@bind it Slider(1:length(tgrid), show_value=true) + +# ╔═╡ 3ea41115-2623-4b55-b5a0-9b3a282b6668 +heatmap(pa[:snaps,1][it], c=:seismic,aspect_ratio=:equal) + # ╔═╡ fd11d02d-80cd-4373-a8d1-e0d50d93dcf8 md""" Lets look at the timer objects. @@ -139,7 +164,16 @@ t1 t2 # ╔═╡ 4612405a-530b-4eab-8b84-40dc907929e4 -plot(plot(data, 99.9), plot(data_box, 99.9), size=(800, 500)) +plot(plot(data, 10), plot(data_box, 10), size=(800, 500)) + +# ╔═╡ c3b9e67c-ec42-422f-bb65-f12269b659d7 +plot(data, 90) + +# ╔═╡ 994331fc-95b7-4417-89b2-bfab6f9b912e +data[:vz] |> maximum + +# ╔═╡ ca5301e3-6740-4f55-a488-121d50e72358 +data[:vx] |> maximum # ╔═╡ ec7985c6-7072-4ac4-9caa-74a815193750 @test data ≠ data_box @@ -160,14 +194,20 @@ plot(plot(data, 99.9), plot(data_box, 99.9), size=(800, 500)) # ╠═10be735e-4012-4fa8-b538-9f9a08694458 # ╠═7fa1ddc5-fe10-44e2-b764-533be1840ab6 # ╟─cec954f6-c278-4003-99e4-6c995f74b606 +# ╠═14713ee7-6753-46a7-98b3-ed0b8ef48ae4 # ╠═c7d96dda-17f5-4dbb-98bf-4c29fd784574 # ╠═9a3360c0-72c7-420f-a7c0-5565a1383f07 # ╠═df4a4072-a91a-471a-88d9-b74ed58e7823 # ╠═07da9113-a518-4925-8978-8407dc60c605 # ╟─7ac4f51d-73b4-4f9c-88a4-a48911742282 # ╠═033cd313-4e87-4364-a1f8-136e0e06bb65 +# ╠═3ea41115-2623-4b55-b5a0-9b3a282b6668 +# ╠═e9c4cb8b-20bf-43c3-8da5-9ecb2e0ddfd2 # ╠═fd11d02d-80cd-4373-a8d1-e0d50d93dcf8 # ╠═368c2bbd-94fb-4826-96eb-8fa37a7e35a3 # ╠═7236dc4a-fe04-492d-9532-170ac7c22fb8 # ╠═4612405a-530b-4eab-8b84-40dc907929e4 +# ╠═c3b9e67c-ec42-422f-bb65-f12269b659d7 +# ╠═994331fc-95b7-4417-89b2-bfab6f9b912e +# ╠═ca5301e3-6740-4f55-a488-121d50e72358 # ╠═ec7985c6-7072-4ac4-9caa-74a815193750 diff --git a/notebooks/sample_fwi.jl b/notebooks/sample_fwi.jl new file mode 100644 index 0000000..5187b29 --- /dev/null +++ b/notebooks/sample_fwi.jl @@ -0,0 +1,502 @@ +### A Pluto.jl notebook ### +# v0.19.21 + +using Markdown +using InteractiveUtils + +# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error). +macro bind(def, element) + quote + local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end + local el = $(esc(element)) + global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) + el + end +end + +# ╔═╡ 474f6970-9784-45fa-9d9b-bb1abbf71929 +using Pkg; Pkg.activate("GeoPhyInv") + +# ╔═╡ 7d9b2219-2074-4cab-8037-9391ee0406fd +Pkg.add("SegyIO"); using SegyIO + +# ╔═╡ 75ad24f7-8b57-401b-8170-d3e35f3c960f +begin + using PlutoLinks: @revise + using PlutoUI, PlutoTest, Plots, CUDA +end + +# ╔═╡ f86410b5-87f7-4761-91a8-856406d00234 +using GeoPhyInv + +# ╔═╡ bd68d9b1-633b-4c9b-b759-791581714306 +using Statistics, LossFunctions, LinearAlgebra, MLUtils, FiniteDifferences, SparseArrays + +# ╔═╡ 903d3dc9-b4d6-44aa-97db-5e36ccbb9f98 +using DSP + +# ╔═╡ 7c9cba9d-7e8f-446f-a664-2e8c988f2095 +using CSV, DataFrames + +# ╔═╡ bbfa8c8e-1ef5-459b-b5f5-edb15ff38fa8 +TableOfContents() + +# ╔═╡ 893cd92d-3d0d-4fb2-b465-f61a237fc154 +# begin +# Pkg.activate(Base.current_project()) +# Pkg.instantiate() +# end + +# ╔═╡ 082b3d0f-bb73-4691-8442-d6dec1310cfd +# GeoPhyInv.set_preferences(ndims=2,use_gpu=true,datatype="Float32",order=2) + +# ╔═╡ 16acf176-4ffa-40f0-a4ce-8d7d8483b5cb + + +# ╔═╡ c60027b8-ff8e-4671-aef3-36e136096b4d +begin + # pa_mod = SeisForwExpt(FdtdAcoustic{FullWave}(:forward, 2), Homogeneous()) + d= 0.1 + mgrid1 = [range(000.0, stop =50.0, step=d), range(0.0, stop = 50, step=d)] + newm=AcousticMedium(mgrid1) + newm1=AcousticMedium(mgrid1) + # pa_mod= + # m1=GeoPhyInv.get_modelvector(pa_mod, [:invK]); + # @time update!(pa_mod, m1, [:invK]) + # m2=GeoPhyInv.get_modelvector(pa_mod, [:invK]) + # @test m1 ≈ m2 +end + +# ╔═╡ 46bdb8a0-350e-41e5-8637-a784b17a81e1 +begin + newm.vp=GeoPhyInv.vp(Float32.(100.0 .* ones(size(newm.vp.m)))) + # medium.vs=GeoPhyInv.vs(Float32.(vs1)) + newm.rho=GeoPhyInv.rho(Float32.(1900.0 .* ones(size(newm.vp.m)))) +end + +# ╔═╡ 33f69a08-639c-4d3f-930d-fe6743d01a3c +begin + newm1.vp=GeoPhyInv.vp(Float32.(100.0 .* ones(size(newm.vp.m)))) + # medium.vs=GeoPhyInv.vs(Float32.(vs1)) + newm1.rho=GeoPhyInv.rho(Float32.(1900.0 .* ones(size(newm.vp.m)))) +end + +# ╔═╡ 09772196-4227-4a7b-80e2-072f3a31231c +GeoPhyInv.update!(newm, [:vp], rectangle = [[12,12], [14,48]], perc =-60) + +# ╔═╡ 8d458409-0d0b-4e41-a0f2-21338377a9f3 +begin + nr=24 + function get_ageom(medium) + ageom=AGeom(medium.grid, :surf, SSrcs(1), Recs(nr)) # surface acquisition + GeoPhyInv.update!(ageom, Recs(nr), [0,10],[0,33]) + # srcs_x =rand(Uniform(10,33)) + srcs_x=10 + # srcs_y = 5 + # srcs_z = rand(Uniform(100,800)) + srcs_z = 0 + + # srcs_x=-29.709914448408938 + # srcs_y=-80.23644969547148 + # println([srcs_x, srcs_y]) + # srcs_x = 0.0 + # srcs_y = 0.0 + # srcs_z=0 + GeoPhyInv.update!(ageom,SSrcs(1), [srcs_z,srcs_x], [srcs_z,srcs_x]) + end + ageom=get_ageom(newm) + ageom1=get_ageom(newm1) + println(ageom) +end + +# ╔═╡ 2fb4a69e-d0b9-4516-a14b-f96a9a37904d +begin + heatmap(mgrid1[2],mgrid1[1],newm1[:vp].m,seriescolor=cgrad(:Blues),yflip=true) + # ,clim=(200,400)) + Plots.scatter!(ageom1,SSrcs(1)) + Plots.scatter!(ageom1,Recs(24)) +end + +# ╔═╡ 4662c24a-24e1-4080-abab-0bc42ff7e39c + + +# ╔═╡ 0e932217-beef-4877-b803-70af59bac5dd +begin + dt=5e-04 + itchunk=Int64(floor(2/dt)) + dt1=5e-4 + rf1=(1/dt)/(1/dt1) +end + +# ╔═╡ dae779e0-3cce-4e72-9ddf-64fdc55858a4 +begin + tgrid=range(dt,stop=2,step=dt) + wav = ricker(25.0, tgrid, tpeak =.5); + wav1 = ricker(4.0, tgrid, tpeak =.5); + # dt=step(tgrid) + # println(dt) + function get_srcwav(wav,i=10.0) + # 8 wavelengths along the diagonal of the medium, and total time of 1 diagonal + + # nts_perc=i + # a=randn(floor(Int,length(tgrid)*nts_perc*0.01)) + # wavv=DSP.conv(a,wav)[1:length(tgrid)] + wavv=wav + # wavv=randobs(read_data(selected_files[randobs(1:size(selected_files,1))]))[1:size(tgrid,1),1] + global freqmax=GeoPhyInv.Utils.findfreq(wavv,tgrid,attrib=:max) + global freqpeak=GeoPhyInv.Utils.findfreq(wavv,tgrid,attrib=:peak) + global freqmin=GeoPhyInv.Utils.findfreq(wavv,tgrid,attrib=:min) + tmax=maximum(tgrid); nt=length(tgrid); + # println("nt:\t",length(tgrid)) + # println("freqpeak:\t",freqpeak) + + + # intialize SrcWavelets + # Srcs(tgrid, ageom, [:p, :vx]) + srcwav = Srcs(tgrid, ageom, [:vz]) + GeoPhyInv.update!(srcwav, [:vz], wavv) + return srcwav + end + srcwav=get_srcwav(wav) + srcwav1=get_srcwav(wav1,2.0) + plot(srcwav[1].d) +end + +# ╔═╡ b0e1670a-2e83-4e3b-86b0-83425169cc3e +md""" +## Load Data +#### Scan a given folder +We provide two options, where the second option overrides the first. + +We search of relevant files in a given input folder. + +$(@bind mainfldrname TextField((50,2), default=".", placeholder="Enter the relative folder path to scan")) + +Select the file extension to scan +$(@bind file_ext Select(["csv", "h5", "dat","sg2","sgy"], default="sgy")) + +#### Forcefully load a particular file +Alternatively, we manually enter the path of the file to read. +$(@bind forced_file FilePicker()) +""" + +# ╔═╡ 18e4bd74-bed3-4f58-8315-f6251b67d943 +@bind selected_folder Select(readdir(mainfldrname),default="Synthetic_data") + +# ╔═╡ 9a44176d-3bb7-4516-95b7-278e03fa3d4a + function filesUI() + if(forced_file===nothing) + filtered_files=filter(readdir(joinpath(mainfldrname,selected_folder), join=true)) do f + (lowercase(splitext(f)[end]) == string(".",file_ext)) + end + return md""" + Depending on the choices made above, select files to preprocess and analyze. By default, 15 files will be randomly selected. + $(@bind selected_files MultiCheckBox(filtered_files, select_all=true, default=sort(first(filtered_files,15))))) + """ + else + return [forced_file] + end +end + +# ╔═╡ 074aad69-bca9-4866-9579-ddec53053eb5 +filesUI() + +# ╔═╡ 2e22d0a7-b0ae-4a22-8aea-439a102b473f +function read_file1(fname, file_ext) +if(file_ext == "sgy") + return Float32.(segy_read(fname).data[:,1:24]) + end +if(file_ext == "h5") + A=Array(first(first(h5open(fname))))[:,:,1] + return A +end + if(file_ext=="csv") + b=CSV.read(fname,DataFrame) + a=[] + for i in 1:24 + if i==1 + a=b[:,i] + else + a=cat(a,b[:,i],dims=2) + end + end + return a + end + # Similarly, add hdf5 and csv support later... +end + +# ╔═╡ a192f25d-0305-4f82-809f-13f16c301084 +Dobs=data=broadcast(selected_files) do fname + read_file1(fname, file_ext) +end + +# ╔═╡ 6bc5b99b-0f30-46c4-94b1-ce7c63cabdab +heatmap(Dobs[1][:, 2:24], clim=(-0.1, 0.1), yflip=true, ylim=(0, 800), c=:seismic) + +# ╔═╡ cc05fa9a-c2d0-445a-ba4a-11530f07f5b3 +heatmap(Dobs[5]) + +# ╔═╡ de9a4882-cba9-4707-a900-3d82b9a2faa1 +begin +# # typically you will run this cell once +# # typically you will run this cell once + tsnaps=tgrid[1:div(length(tgrid),20):end] # store 20 snapshots + rfields=[:vz] # add as many fields are possible here to record +# pa_true1 = SeisForwExpt( +# #FdtdAcou(), +# FdtdAcoustic{FullWave}(:forward, 2), +# snaps_field=:vz, # comment these if not testing +# tsnaps=tsnaps, # comment these if not testing +# # # pml_edges=[:xmin, :xmax, :zmax, :ymin, :ymax], # if you want surface waves to be modelled +# pml_faces=[:zmax, :xmax, :xmin], # absorbing BC on edges +# ageom = ageom, + +# srcwav = srcwav, +# medium = newm, +# # rfields=rfields, +# rfields=rfields, + +# # sflags = 1, +# # rflags = 1, +# tgrid = tgrid, +# stressfree_faces=[:zmin], # location of the stree-free boundary +# verbose = true, +# ); +end + +# ╔═╡ 1f0eea14-c416-4bdf-83d8-908fc0d8649d +# begin +# # pa_true = SeisForwExpt(FdtdAcoustic{FullWave}(:forward, 1), RandScatterer()) +# update!(pa_true1) +# dobs = deepcopy(pa_true1.c.data[1]) +# end + +# ╔═╡ 2a737f8b-92d0-4042-96eb-6b85ea4030c7 +# @time update!(m1, pa_mod, [:invK]) + +# ╔═╡ 000a9016-2b39-428b-a635-3accd26101a1 +begin +# heatmap(pa_true1[:data][:vz],yflip=true, title="vz", c=:seismic) +# p2=heatmap(pa[:data][:vy], title="vy", c=:seismic) +# p3=heatmap(pa_mod[:data][:vx],yflip=true, title="vx", c=:seismic) +# plot(p1,p3, size=(1000,400), layout=(1,2)) +end + +# ╔═╡ a715feb4-8270-40a4-8796-da1feb6e972e + pa_mod = SeisForwExpt( + #FdtdAcou(), + FdtdAcoustic{FullWave}(:forward, 2), + # snaps_field=:vz, # comment these if not testing + # tsnaps=tsnaps, # comment these if not testing + # # pml_edges=[:xmin, :xmax, :zmax, :ymin, :ymax], # if you want surface waves to be modelled + pml_faces=[:zmax, :xmax, :xmin], # absorbing BC on edges + ageom = ageom1, + srcwav = srcwav, + medium = newm1, + # rfields=rfields, + rfields=rfields, + + # sflags = 1, + # rflags = 1, + tgrid = tgrid, + stressfree_faces=[:zmin], # location of the stree-free boundary + verbose = true, + ); + +# ╔═╡ a35aefc6-4de0-4886-8366-acda8b92e22e +update!(pa_mod) + +# ╔═╡ 839dd5fb-e555-4872-b6f4-342f501a4e78 +begin + pa_true=deepcopy(pa_mod) + update!(pa_true, newm) + update!(pa_true) +end + +# ╔═╡ c1751f48-7c26-4bac-8a47-77bf9c4b3832 +plot(pa_true.c.data[1][1], 99) + +# ╔═╡ a2618c5d-5c6a-4a73-b1a6-cbafb794f224 +# pa_true.c.data[1][1][:vz].=data[1] + +# ╔═╡ 0cf61f2f-435b-4fea-b6f7-e7458cb06b89 +# update!(pa_true) + +# ╔═╡ 78cd9487-b705-4a64-b53b-2caa13946003 +begin +heatmap(pa_mod[:data][:vz],yflip=true, title="vz", c=:seismic) +# p2=heatmap(pa[:data][:vy], title="vy", c=:seismic) +# p3=heatmap(pa_mod[:data][:vx],yflip=true, title="vx", c=:seismic) +# plot(p1,p3, size=(1000,400), layout=(1,2)) +end + +# ╔═╡ df0ff9a7-d7ab-4d9d-a043-d8359d75b1da +plot(pa_mod[:data]) + +# ╔═╡ f49439df-3db3-48ae-bf9d-1387e3c14c8f +N = 10 + +# ╔═╡ 6d47ae1d-2554-4d42-af6c-02ce952bf14e +@bind mpara Select([:invK, :rho]) + +# ╔═╡ 5582110f-327c-419a-8759-ca11844c0c96 +mpara + +# ╔═╡ 51030579-4d25-4e03-ae49-38071dd58a40 +dobs = deepcopy(pa_mod.c.data[1]) + +# ╔═╡ 9d566f59-fba0-4bd5-9945-49566bf74645 +dobs[1][:vz].=Dobs[1] + +# ╔═╡ 26c9f7a1-f18b-4fe3-a607-ad1f6a904bb8 +pa_inv, pa_src = SeisInvExpt(pa_mod, dobs, [range(0,50,length=N), range(0,50,length=N)], [mpara]) + +# ╔═╡ 35d6fbbc-95a9-443b-924a-8f1946e6bcb8 +heatmap(Dobs[1]) + +# ╔═╡ 0e5eb59a-aad6-4c84-a0e4-c8558964719f +# update!(pa_inv, pa_src) + +# ╔═╡ e57e1a51-d69b-4459-9ae9-8869782a0721 +# pa_inv + +# ╔═╡ 1b102fdb-3e82-4b11-b100-fbe1b0588543 +m2 = GeoPhyInv.get_modelvector(pa_inv) + +# ╔═╡ 17ea5b1d-4296-47e0-a361-aaffb3c152f8 +iszero(m2) + +# ╔═╡ 26b174f6-3fbc-4f02-9279-bd83908b1aa7 +LV(m)=GeoPhyInv.lossvalue(m,pa_inv) + +# ╔═╡ d015a19a-a5bd-4660-ae22-dde2ae179313 +Grad(g,m)=GeoPhyInv.gradient!(g, m, pa_inv) + +# ╔═╡ ccbe14ab-5bc6-43ab-8e7b-14076cf37340 +gtest=similar(m2) + +# ╔═╡ cea31a54-6e9c-4505-ac40-8e850f98f9f1 +GeoPhyInv.gradient!(gtest,m2,pa_inv) + +# ╔═╡ 2d9177b5-c10d-4027-8b88-574863609031 +# Grad(gtest,m2) + +# ╔═╡ f470fc4b-4506-4039-a21f-2d4e9ae67e2a +LV(m2) + +# ╔═╡ 406f05a3-4691-4405-a272-579d44aa3ff7 +pa_inv + +# ╔═╡ 7a956543-5a19-4912-958c-521b3ecc7390 +length(pa_mod.c.exmedium.grid[1]) + +# ╔═╡ 2d0eda58-022e-4836-bde8-e4af663d20fb +pa_inv.P +# gtest + +# ╔═╡ 6ec9f2b1-be56-411a-8f3c-518e80d9c078 +heatmap(reshape(Array(pa_inv.mfull),length(pa_mod.c.exmedium.grid[1]),length(pa_mod.c.exmedium.grid[2]))) + +# ╔═╡ 654bf0b4-21d4-419c-bf19-c7909e6b2134 +begin + heatmap(mgrid1[2],mgrid1[1],newm[:vp].m,seriescolor=cgrad(:Blues),yflip=true) + # ,clim=(200,400)) + Plots.scatter!(ageom,SSrcs()) + Plots.scatter!(ageom,Recs(24)) +end + +# ╔═╡ 52c54384-458b-410f-b393-10d475d4f13d +GeoPhyInv.get_modelvector(pa_inv) + +# ╔═╡ 12a6565b-3edd-4e8b-8895-eb49366942ec +heatmap(pa_mod.c.exmedium.grid[2],pa_mod.c.exmedium.grid[1],reshape(Array(pa_inv.P'*GeoPhyInv.get_modelvector(pa_inv)),length(pa_mod.c.exmedium.grid[1]),length(pa_mod.c.exmedium.grid[2])),yflip=true,c=:seismic) + +# ╔═╡ 30598e33-a52e-4c14-895e-7ceb45865573 +heatmap(pa_mod.c.exmedium.grid[2],pa_mod.c.exmedium.grid[1],reshape(Array(pa_inv.P'*gtest),length(pa_mod.c.exmedium.grid[1]),length(pa_mod.c.exmedium.grid[2])),yflip=true,c=:seismic) + +# ╔═╡ db243c16-e57e-4ea4-a291-eab647188023 +heatmap(pa_mod.c.exmedium.grid[2],pa_mod.c.exmedium.grid[1],Array(deepcopy(pa_mod.c.gradients[:invK])),yflip=true,c=:seismic) + +# ╔═╡ 28d9ed78-fd78-46f6-8fd1-13ebfedf1561 +heatmap(pa_true.c.exmedium.grid[2],pa_true.c.exmedium.grid[1],Array(deepcopy(pa_true.c.gradients[:invK])),clim=(-1e-12,1e-12),yflip=true,c=:seismic) + +# ╔═╡ a11ea26a-ab85-468b-809e-34144eb6578f +# using Optim + +# ╔═╡ 6180ba9a-5aee-4e42-918c-5b6bb7728655 +# using Ipopt,JuMP + +# ╔═╡ 43fb3211-9ffc-4ea5-a459-539cf55ba009 +# pa_inv = SeisInvExpt(pa_mod, dobs, [N, N], [mpara]) + +# ╔═╡ Cell order: +# ╠═bbfa8c8e-1ef5-459b-b5f5-edb15ff38fa8 +# ╠═474f6970-9784-45fa-9d9b-bb1abbf71929 +# ╠═75ad24f7-8b57-401b-8170-d3e35f3c960f +# ╠═893cd92d-3d0d-4fb2-b465-f61a237fc154 +# ╠═f86410b5-87f7-4761-91a8-856406d00234 +# ╠═7d9b2219-2074-4cab-8037-9391ee0406fd +# ╠═082b3d0f-bb73-4691-8442-d6dec1310cfd +# ╠═bd68d9b1-633b-4c9b-b759-791581714306 +# ╠═16acf176-4ffa-40f0-a4ce-8d7d8483b5cb +# ╠═c60027b8-ff8e-4671-aef3-36e136096b4d +# ╠═46bdb8a0-350e-41e5-8637-a784b17a81e1 +# ╠═33f69a08-639c-4d3f-930d-fe6743d01a3c +# ╠═09772196-4227-4a7b-80e2-072f3a31231c +# ╠═8d458409-0d0b-4e41-a0f2-21338377a9f3 +# ╠═2fb4a69e-d0b9-4516-a14b-f96a9a37904d +# ╠═4662c24a-24e1-4080-abab-0bc42ff7e39c +# ╠═0e932217-beef-4877-b803-70af59bac5dd +# ╠═903d3dc9-b4d6-44aa-97db-5e36ccbb9f98 +# ╠═dae779e0-3cce-4e72-9ddf-64fdc55858a4 +# ╠═b0e1670a-2e83-4e3b-86b0-83425169cc3e +# ╠═18e4bd74-bed3-4f58-8315-f6251b67d943 +# ╠═074aad69-bca9-4866-9579-ddec53053eb5 +# ╠═9a44176d-3bb7-4516-95b7-278e03fa3d4a +# ╠═2e22d0a7-b0ae-4a22-8aea-439a102b473f +# ╠═7c9cba9d-7e8f-446f-a664-2e8c988f2095 +# ╠═a192f25d-0305-4f82-809f-13f16c301084 +# ╠═6bc5b99b-0f30-46c4-94b1-ce7c63cabdab +# ╠═cc05fa9a-c2d0-445a-ba4a-11530f07f5b3 +# ╠═de9a4882-cba9-4707-a900-3d82b9a2faa1 +# ╠═1f0eea14-c416-4bdf-83d8-908fc0d8649d +# ╠═2a737f8b-92d0-4042-96eb-6b85ea4030c7 +# ╠═000a9016-2b39-428b-a635-3accd26101a1 +# ╠═a715feb4-8270-40a4-8796-da1feb6e972e +# ╠═a35aefc6-4de0-4886-8366-acda8b92e22e +# ╠═839dd5fb-e555-4872-b6f4-342f501a4e78 +# ╠═c1751f48-7c26-4bac-8a47-77bf9c4b3832 +# ╠═a2618c5d-5c6a-4a73-b1a6-cbafb794f224 +# ╠═0cf61f2f-435b-4fea-b6f7-e7458cb06b89 +# ╠═78cd9487-b705-4a64-b53b-2caa13946003 +# ╠═df0ff9a7-d7ab-4d9d-a043-d8359d75b1da +# ╠═f49439df-3db3-48ae-bf9d-1387e3c14c8f +# ╠═6d47ae1d-2554-4d42-af6c-02ce952bf14e +# ╠═5582110f-327c-419a-8759-ca11844c0c96 +# ╠═51030579-4d25-4e03-ae49-38071dd58a40 +# ╠═9d566f59-fba0-4bd5-9945-49566bf74645 +# ╠═26c9f7a1-f18b-4fe3-a607-ad1f6a904bb8 +# ╠═35d6fbbc-95a9-443b-924a-8f1946e6bcb8 +# ╠═0e5eb59a-aad6-4c84-a0e4-c8558964719f +# ╠═e57e1a51-d69b-4459-9ae9-8869782a0721 +# ╠═1b102fdb-3e82-4b11-b100-fbe1b0588543 +# ╠═17ea5b1d-4296-47e0-a361-aaffb3c152f8 +# ╠═26b174f6-3fbc-4f02-9279-bd83908b1aa7 +# ╠═d015a19a-a5bd-4660-ae22-dde2ae179313 +# ╠═ccbe14ab-5bc6-43ab-8e7b-14076cf37340 +# ╠═cea31a54-6e9c-4505-ac40-8e850f98f9f1 +# ╠═2d9177b5-c10d-4027-8b88-574863609031 +# ╠═f470fc4b-4506-4039-a21f-2d4e9ae67e2a +# ╠═406f05a3-4691-4405-a272-579d44aa3ff7 +# ╠═7a956543-5a19-4912-958c-521b3ecc7390 +# ╠═2d0eda58-022e-4836-bde8-e4af663d20fb +# ╠═6ec9f2b1-be56-411a-8f3c-518e80d9c078 +# ╠═654bf0b4-21d4-419c-bf19-c7909e6b2134 +# ╠═52c54384-458b-410f-b393-10d475d4f13d +# ╠═12a6565b-3edd-4e8b-8895-eb49366942ec +# ╠═30598e33-a52e-4c14-895e-7ceb45865573 +# ╠═db243c16-e57e-4ea4-a291-eab647188023 +# ╠═28d9ed78-fd78-46f6-8fd1-13ebfedf1561 +# ╠═a11ea26a-ab85-468b-809e-34144eb6578f +# ╠═6180ba9a-5aee-4e42-918c-5b6bb7728655 +# ╠═43fb3211-9ffc-4ea5-a459-539cf55ba009 diff --git a/src/GeoPhyInv.jl b/src/GeoPhyInv.jl index acbeb8b..ccc5338 100644 --- a/src/GeoPhyInv.jl +++ b/src/GeoPhyInv.jl @@ -37,7 +37,6 @@ using StatsBase using InteractiveUtils using LaTeXStrings using RecipesBase -using ColorSchemes using FFTW using HDF5 using SpecialFunctions @@ -88,8 +87,8 @@ const _fd_order = @load_preference("order", 2) const _fd_ndims = @load_preference("ndims", 2) const _fd_datatype = eval(Symbol(@load_preference("datatype", "Float32"))) const _fd_use_gpu = @load_preference("use_gpu", false) -const _fd_npml = 20 + (_fd_order - 1) -const _fd_npextend = 20 + (_fd_order - 1) # determines exmedium +const _fd_npml = 40 + (_fd_order - 1) +const _fd_npextend = 40 + (_fd_order - 1) # determines exmedium const _fd_nbound = 3 # number of points to store the boundary values diff --git a/src/fdtd/advance_elastic.jl b/src/fdtd/advance_elastic.jl index c64a7d8..72be293 100644 --- a/src/fdtd/advance_elastic.jl +++ b/src/fdtd/advance_elastic.jl @@ -165,22 +165,22 @@ end @all(tauxx) = @all(tauxx) - - (@all(dtM) * @all(dvxdx)) + (@all(dtlambda) * (@all(dvydy) + @all(dvzdz))) + (@all(dtM) * @all(dvxdx)) - (@all(dtlambda) * (@all(dvydy) + @all(dvzdz))) @all(tauyy) = @all(tauyy) - - (@all(dtM) * @all(dvydy)) + (@all(dtlambda) * (@all(dvxdx) + @all(dvzdz))) + (@all(dtM) * @all(dvydy)) - (@all(dtlambda) * (@all(dvxdx) + @all(dvzdz))) @all(tauzz) = @all(tauzz) - - (@all(dtM) * @all(dvzdz)) + (@all(dtlambda) * (@all(dvydy) + @all(dvxdx))) + (@all(dtM) * @all(dvzdz)) - (@all(dtlambda) * (@all(dvydy) + @all(dvxdx))) return end @parallel function compute_stressii!(tauxx::Data.Array{2}, tauzz, dvxdx, dvzdz, dtM, dtlambda) @all(tauxx) = - @all(tauxx) - (@all(dtM) * @all(dvxdx)) + (@all(dtlambda) * (@all(dvzdz))) + @all(tauxx) - (@all(dtM) * @all(dvxdx)) - (@all(dtlambda) * (@all(dvzdz))) @all(tauzz) = - @all(tauzz) - (@all(dtM) * @all(dvzdz)) + (@all(dtlambda) * (@all(dvxdx))) + @all(tauzz) - (@all(dtM) * @all(dvzdz)) - (@all(dtlambda) * (@all(dvxdx))) return end diff --git a/src/fdtd/diff2D.jl b/src/fdtd/diff2D.jl index d697a56..61056a9 100644 --- a/src/fdtd/diff2D.jl +++ b/src/fdtd/diff2D.jl @@ -1,6 +1,6 @@ # borrowed 2-point stencil macros from ParallelStencil.FiniteDifferences; extended to higher-order differences -import ..ParallelStencil: INDICES, WITHIN_DOC +import ..ParallelStencil: INDICES, WITHIN_DOC, @expandargs iz, ix = INDICES[1], INDICES[2] @static if (_fd_order == 2) izi, ixi = :($iz + 1), :($ix + 1) @@ -15,6 +15,7 @@ end @static if (_fd_order == 2) macro within(macroname::String, A::Symbol) + @expandargs(A) if macroname == "@all" esc(:($iz <= size($A, 1) && $ix <= size($A, 2))) elseif macroname == "@inn" @@ -27,21 +28,26 @@ end end macro d_za(A::Symbol) + @expandargs(A) esc(:($A[$iz+1, $ix] - $A[$iz, $ix])) end macro d_xa(A::Symbol) + @expandargs(A) esc(:($A[$iz, $ix+1] - $A[$iz, $ix])) end macro d_zi(A::Symbol) + @expandargs(A) esc(:($A[$iz+1, $ixi] - $A[$iz, $ixi])) end macro d_xi(A::Symbol) + @expandargs(A) esc(:($A[$izi, $ix+1] - $A[$izi, $ix])) end # dummy y elseif (_fd_order == 4) macro within(macroname::String, A::Symbol) + @expandargs(A) if macroname == "@all" esc(:($iz <= size($A, 1) && $ix <= size($A, 2))) elseif macroname == "@inn" @@ -54,6 +60,7 @@ elseif (_fd_order == 4) end macro d_za(A::Symbol) + @expandargs(A) esc( :( $A[$iz+2, $ix] * 27.0 - $A[$iz+1, $ix] * 27.0 + $A[$iz, $ix] - @@ -62,6 +69,7 @@ elseif (_fd_order == 4) ) end macro d_xa(A::Symbol) + @expandargs(A) esc( :( $A[$iz, $ix+2] * 27.0 - $A[$iz, $ix+1] * 27.0 + $A[$iz, $ix] - @@ -70,6 +78,7 @@ elseif (_fd_order == 4) ) end macro d_zi(A::Symbol) + @expandargs(A) esc( :( $A[$iz+2, $ixi] * 27.0 - $A[$iz+1, $ixi] * 27.0 + $A[$iz, $ixi] - @@ -78,6 +87,7 @@ elseif (_fd_order == 4) ) end macro d_xi(A::Symbol) + @expandargs(A) esc( :( $A[$izi, $ix+2] * 27.0 - $A[$izi, $ix+1] * 27.0 + $A[$izi, $ix] - @@ -87,6 +97,7 @@ elseif (_fd_order == 4) end elseif (_fd_order == 6) macro within(macroname::String, A::Symbol) + @expandargs(A) if macroname == "@all" esc(:($iz <= size($A, 1) && $ix <= size($A, 2))) elseif macroname == "@inn" @@ -99,6 +110,7 @@ elseif (_fd_order == 6) end macro d_za(A::Symbol) + @expandargs(A) esc( :( $A[$iz+3, $ix] * 2250.0 - $A[$iz+2, $ix] * 2250.0 + $A[$iz+1, $ix] * 125.0 - @@ -107,6 +119,7 @@ elseif (_fd_order == 6) ) end macro d_xa(A::Symbol) + @expandargs(A) esc( :( $A[$iz, $ix+3] * 2250.0 - $A[$iz, $ix+2] * 2250.0 + $A[$iz, $ix+1] * 125.0 - @@ -115,6 +128,7 @@ elseif (_fd_order == 6) ) end macro d_zi(A::Symbol) + @expandargs(A) esc( :( $A[$iz+3, $ixi] * 2250.0 - $A[$iz+2, $ixi] * 2250.0 + @@ -124,6 +138,7 @@ elseif (_fd_order == 6) ) end macro d_xi(A::Symbol) + @expandargs(A) esc( :( $A[$izi, $ix+3] * 2250.0 - $A[$izi, $ix+2] * 2250.0 + @@ -133,6 +148,7 @@ elseif (_fd_order == 6) ) end elseif (_fd_order == 8) + @expandargs(A) macro within(macroname::String, A::Symbol) if macroname == "@all" esc(:($iz <= size($A, 1) && $ix <= size($A, 2))) @@ -146,6 +162,7 @@ elseif (_fd_order == 8) end macro d_za(A::Symbol) + @expandargs(A) esc( :( $A[$iz+4, $ix] * 128625.0 - $A[$iz+3, $ix] * 128625.0 + @@ -156,6 +173,7 @@ elseif (_fd_order == 8) ) end macro d_xa(A::Symbol) + @expandargs(A) esc( :( $A[$iz, $ix+4] * 128625.0 - $A[$iz, $ix+3] * 128625.0 + @@ -166,6 +184,7 @@ elseif (_fd_order == 8) ) end macro d_zi(A::Symbol) + @expandargs(A) esc( :( $A[$iz+4, $ixi] * 128625.0 - $A[$iz+3, $ixi] * 128625.0 + @@ -176,6 +195,7 @@ elseif (_fd_order == 8) ) end macro d_xi(A::Symbol) + @expandargs(A) esc( :( $A[$izi, $ix+4] * 128625.0 - $A[$izi, $ix+3] * 128625.0 + @@ -192,18 +212,23 @@ elseif (_fd_order == 8) end macro all(A::Symbol) + @expandargs(A) esc(:($A[$iz, $ix])) end macro inn(A::Symbol) + @expandargs(A) esc(:($A[$izi, $ixi])) end macro av(A::Symbol) + @expandargs(A) esc(:(($A[$iz, $ix] + $A[$iz+1, $ix] + $A[$iz, $ix+1] + $A[$iz+1, $ix+1]) * 0.25)) end macro av_zi(A::Symbol) + @expandargs(A) esc(:(($A[$iz, $ixi] + $A[$iz+1, $ixi]) * 0.5)) end macro av_xi(A::Symbol) + @expandargs(A) esc(:(($A[$izi, $ix] + $A[$izi, $ix+1]) * 0.5)) end diff --git a/src/fdtd/diff3D.jl b/src/fdtd/diff3D.jl index 3394c3e..eaf342c 100644 --- a/src/fdtd/diff3D.jl +++ b/src/fdtd/diff3D.jl @@ -1,6 +1,6 @@ # borrowed 2-point stencil macros from ParallelStencil.FiniteDifferences; extended to higher-order differences -import ..ParallelStencil: INDICES, WITHIN_DOC +import ..ParallelStencil: INDICES, WITHIN_DOC, @expandargs iz, iy, ix = INDICES[1], INDICES[2], INDICES[3] @static if (_fd_order == 2) izi, iyi, ixi = :($iz + 1), :($iy + 1), :($ix + 1) @@ -15,6 +15,7 @@ end @static if (_fd_order == 2) macro within(macroname::String, A::Symbol) + @expandargs(A) if macroname == "@all" esc(:($iz <= size($A, 1) && $iy <= size($A, 2) && $ix <= size($A, 3))) elseif macroname == "@inn" @@ -33,26 +34,33 @@ end end macro d_za(A::Symbol) + @expandargs(A) esc(:($A[$iz+1, $iy, $ix] - $A[$iz, $iy, $ix])) end macro d_ya(A::Symbol) + @expandargs(A) esc(:($A[$iz, $iy+1, $ix] - $A[$iz, $iy, $ix])) end macro d_xa(A::Symbol) + @expandargs(A) esc(:($A[$iz, $iy, $ix+1] - $A[$iz, $iy, $ix])) end macro d_zi(A::Symbol) + @expandargs(A) esc(:($A[$iz+1, $iyi, $ixi] - $A[$iz, $iyi, $ixi])) end macro d_yi(A::Symbol) + @expandargs(A) esc(:($A[$izi, $iy+1, $ixi] - $A[$izi, $iy, $ixi])) end macro d_xi(A::Symbol) + @expandargs(A) esc(:($A[$izi, $iyi, $ix+1] - $A[$izi, $iyi, $ix])) end elseif (_fd_order == 4) macro within(macroname::String, A::Symbol) + @expandargs(A) if macroname == "@all" esc(:($iz <= size($A, 1) && $iy <= size($A, 2) && $ix <= size($A, 3))) elseif macroname == "@inn" @@ -71,6 +79,7 @@ elseif (_fd_order == 4) end macro d_za(A::Symbol) + @expandargs(A) esc( :( $A[$iz+2, $iy, $ix] * 27.0 - $A[$iz+1, $iy, $ix] * 27.0 + @@ -79,6 +88,7 @@ elseif (_fd_order == 4) ) end macro d_ya(A::Symbol) + @expandargs(A) esc( :( $A[$iz, $iy+2, $ix] * 27.0 - $A[$iz, $iy+1, $ix] * 27.0 + @@ -87,6 +97,7 @@ elseif (_fd_order == 4) ) end macro d_xa(A::Symbol) + @expandargs(A) esc( :( $A[$iz, $iy, $ix+2] * 27.0 - $A[$iz, $iy, $ix+1] * 27.0 + @@ -95,6 +106,7 @@ elseif (_fd_order == 4) ) end macro d_zi(A::Symbol) + @expandargs(A) esc( :( $A[$iz+2, $iyi, $ixi] * 27.0 - $A[$iz+1, $iyi, $ixi] * 27.0 + @@ -103,6 +115,7 @@ elseif (_fd_order == 4) ) end macro d_yi(A::Symbol) + @expandargs(A) esc( :( $A[$izi, $iy+2, $ixi] * 27.0 - $A[$izi, $iy+1, $ixi] * 27.0 + @@ -111,6 +124,7 @@ elseif (_fd_order == 4) ) end macro d_xi(A::Symbol) + @expandargs(A) esc( :( $A[$izi, $iyi, $ix+2] * 27.0 - $A[$izi, $iyi, $ix+1] * 27.0 + @@ -120,6 +134,7 @@ elseif (_fd_order == 4) end elseif (_fd_order == 6) macro within(macroname::String, A::Symbol) + @expandargs(A) if macroname == "@all" esc(:($iz <= size($A, 1) && $ix <= size($A, 2))) elseif macroname == "@inn" @@ -138,6 +153,7 @@ elseif (_fd_order == 6) end macro d_za(A::Symbol) + @expandargs(A) esc( :( $A[$iz+3, $iy, $ix] * 2250.0 - $A[$iz+2, $iy, $ix] * 2250.0 + @@ -147,6 +163,7 @@ elseif (_fd_order == 6) ) end macro d_ya(A::Symbol) + @expandargs(A) esc( :( $A[$iz, $iy+3, $ix] * 2250.0 - $A[$iz, $iy+2, $ix] * 2250.0 + @@ -156,6 +173,7 @@ elseif (_fd_order == 6) ) end macro d_xa(A::Symbol) + @expandargs(A) esc( :( $A[$iz, $iy, $ix+3] * 2250.0 - $A[$iz, $iy, $ix+2] * 2250.0 + @@ -165,6 +183,7 @@ elseif (_fd_order == 6) ) end macro d_zi(A::Symbol) + @expandargs(A) esc( :( $A[$iz+3, $iyi, $ixi] * 2250.0 - $A[$iz+2, $iyi, $ixi] * 2250.0 + @@ -174,6 +193,7 @@ elseif (_fd_order == 6) ) end macro d_yi(A::Symbol) + @expandargs(A) esc( :( $A[$izi, $iy+3, $ixi] * 2250.0 - $A[$izi, $iy+2, $ixi] * 2250.0 + @@ -183,6 +203,7 @@ elseif (_fd_order == 6) ) end macro d_xi(A::Symbol) + @expandargs(A) esc( :( $A[$izi, $iyi, $ix+3] * 2250.0 - $A[$izi, $iyi, $ix+2] * 2250.0 + @@ -193,6 +214,7 @@ elseif (_fd_order == 6) end elseif (_fd_order == 8) macro within(macroname::String, A::Symbol) + @expandargs(A) if macroname == "@all" esc(:($iz <= size($A, 1) && $ix <= size($A, 2))) elseif macroname == "@inn" @@ -211,6 +233,7 @@ elseif (_fd_order == 8) end macro d_za(A::Symbol) + @expandargs(A) esc( :( $A[$iz+4, $iy, $ix] * 128625.0 - $A[$iz+3, $iy, $ix] * 128625.0 + @@ -221,6 +244,7 @@ elseif (_fd_order == 8) ) end macro d_ya(A::Symbol) + @expandargs(A) esc( :( $A[$iz, $iy+4, $ix] * 128625.0 - $A[$iz, $iy+3, $ix] * 128625.0 + @@ -231,6 +255,7 @@ elseif (_fd_order == 8) ) end macro d_xa(A::Symbol) + @expandargs(A) esc( :( $A[$iz, $iy, $ix+4] * 128625.0 - $A[$iz, $iy, $ix+3] * 128625.0 + @@ -241,6 +266,7 @@ elseif (_fd_order == 8) ) end macro d_zi(A::Symbol) + @expandargs(A) esc( :( $A[$iz+4, $iyi, $ixi] * 128625.0 - $A[$iz+3, $iyi, $ixi] * 128625.0 + @@ -251,6 +277,7 @@ elseif (_fd_order == 8) ) end macro d_yi(A::Symbol) + @expandargs(A) esc( :( $A[$izi, $iy+4, $ixi] * 128625.0 - $A[$izi, $iy+3, $ixi] * 128625.0 + @@ -261,6 +288,7 @@ elseif (_fd_order == 8) ) end macro d_xi(A::Symbol) + @expandargs(A) esc( :( $A[$izi, $iyi, $ix+4] * 128625.0 - $A[$izi, $iyi, $ix+3] * 128625.0 + @@ -273,12 +301,15 @@ elseif (_fd_order == 8) end macro all(A::Symbol) + @expandargs(A) esc(:($A[$iz, $iy, $ix])) end macro inn(A::Symbol) + @expandargs(A) esc(:($A[$izi, $iyi, $ixi])) end macro av(A::Symbol) + @expandargs(A) esc( :( ( @@ -295,15 +326,19 @@ macro av(A::Symbol) ) end macro av_zi(A::Symbol) + @expandargs(A) esc(:(($A[$iz, $iyi, $ixi] + $A[$iz+1, $iyi, $ixi]) * 0.5)) end macro av_yi(A::Symbol) + @expandargs(A) esc(:(($A[$izi, $iy, $ixi] + $A[$izi, $iy+1, $ixi]) * 0.5)) end macro av_xi(A::Symbol) + @expandargs(A) esc(:(($A[$izi, $iyi, $ix] + $A[$izi, $iyi, $ix+1]) * 0.5)) end macro av_yzi(A::Symbol) + @expandargs(A) esc( :( ( @@ -316,6 +351,7 @@ macro av_yzi(A::Symbol) ) end macro av_xzi(A::Symbol) + @expandargs(A) esc( :( ( @@ -328,6 +364,7 @@ macro av_xzi(A::Symbol) ) end macro av_xyi(A::Symbol) + @expandargs(A) esc( :( ( diff --git a/src/fdtd/fdtd.jl b/src/fdtd/fdtd.jl index 2feba56..a5fade6 100644 --- a/src/fdtd/fdtd.jl +++ b/src/fdtd/fdtd.jl @@ -322,6 +322,7 @@ function P_x_worker_x_pw(ipw, sschunks::UnitRange{Int64}, pac::P_common{T,N}) wh [:t, :tp], # current and previous time step ) velocity_buffer = NamedArray([zeros(eval(f)(), T(), n...) for f in velocity_fields], Symbol.(velocity_fields)) + tauii_buffer = Data.Array(zeros(n...)) # dummy (use for viscoelastic modeling later) w2 = NamedArray( @@ -338,7 +339,7 @@ function P_x_worker_x_pw(ipw, sschunks::UnitRange{Int64}, pac::P_common{T,N}) wh ss = [P_x_worker_x_pw_x_ss(ipw, iss, pac) for (issp, iss) in enumerate(sschunks)] - return P_x_worker_x_pw(ss, w1, w2, memory_pml, velocity_buffer) + return P_x_worker_x_pw(ss, w1, w2, memory_pml, velocity_buffer, tauii_buffer) end diff --git a/src/fdtd/propagate.jl b/src/fdtd/propagate.jl index d7f5a67..9b4f50c 100644 --- a/src/fdtd/propagate.jl +++ b/src/fdtd/propagate.jl @@ -199,7 +199,7 @@ o o o o o end @timeit to "add body-force velocity sources" begin - add_source!(it, issp, iss, pac, pap, activepw, src_flags, [:vx, :vy, :vz]) + add_velocity_source!(it, issp, iss, pac, pap, activepw, src_flags, [:vx, :vy, :vz]) end add_born_sources_velocity!(Val(pac.attrib_mod.mode), pap, pac) @@ -220,7 +220,7 @@ o o o o o end @timeit to "add stress sources" begin # only pressure source, for now - add_source!(it, issp, iss, pac, pap, activepw, src_flags, [:p]) + add_stress_source!(it, issp, iss, pac, pap, activepw, src_flags, [:p, :tauxx, :tauyy, :tauzz]) end add_born_sources_stress!(Val(pac.attrib_mod.mode), pap, pac) diff --git a/src/fdtd/source.jl b/src/fdtd/source.jl index 3e1428d..75cffac 100644 --- a/src/fdtd/source.jl +++ b/src/fdtd/source.jl @@ -3,7 +3,7 @@ get_source(w, ::Any, ::Val{0}) = zero(w) -get_source(w, ::Union{vz,vx,vy}, ::Val{1}) = copy(w) +get_source(w, ::Union{vz,vx,vy,p,tauxx,tauyy,tauzz}, ::Val{1}) = copy(w) # -ve values correspond to source sink, used while solving boundary value problem function get_source(w, ::Union{vz,vx,vy}, ::Val{-1}) ww = copy(w) @@ -22,7 +22,6 @@ end Use input srcwav to fill wavelets and returns frequency bounds. """ function fill_wavelets!(ipw, iss, wavelets, srcwav, src_types) - nt = size(wavelets, 1) sfields = names(srcwav[ipw][iss].d)[1] ns = srcwav[ipw][iss][:ns] @@ -58,8 +57,71 @@ function fill_wavelets!(ipw, iss, wavelets, srcwav, src_types) return freqmin, freqmax, Statistics.mean(freqpeaks) 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, pap, activepw, src_flags, sfields=Fields(pac.attrib_mod)) + +@inbounds @fastmath function add_stress_source!(it::Int64, issp::Int64, iss::Int64, pac::T, pap, activepw, src_flags, sfields=Fields(pac.attrib_mod)) where {T<:P_common{<:FdtdAcoustic}} + if (1 ∈ activepw) + ipw = 1 + sfields = intersect(names(pac.srcwav[ipw][iss].d)[1], sfields) + for sfield in sfields + if (src_flags[ipw]) # add only if src_flags is non-zero + pv = pap[ipw].tauii_buffer # pick stress buffer + pv_vec = view(pv, :) + w = pap[ipw].ss[issp].wavelets[it][sfield] + ssprayw = pap[ipw].ss[issp].ssprayw[sfield] + mul!(pv_vec, ssprayw, w) + pw = pap[ipw].w1[:t][:p] + @parallel muladd_tauii!(pw, pv, pac.dmod[:dtK]) + end + end + end +end + +@inbounds @fastmath function add_stress_source!(it::Int64, issp::Int64, iss::Int64, pac::T, pap, activepw, src_flags, sfields=Fields(pac.attrib_mod)) where {T<:P_common{<:FdtdElastic,2}} + if (1 ∈ activepw) + ipw = 1 + sfields = intersect(names(pac.srcwav[ipw][iss].d)[1], sfields) + for sfield in sfields + if (src_flags[ipw]) # add only if src_flags is non-zero + pv = pap[ipw].tauii_buffer # pick stress buffer + pv_vec = view(pv, :) + w = pap[ipw].ss[issp].wavelets[it][sfield] + ssprayw = pap[ipw].ss[issp].ssprayw[sfield] + mul!(pv_vec, ssprayw, w) + pw = pap[ipw].w1[:t][:tauxx] + @parallel muladd_tauii!(pw, pv, pac.dmod[:dtM]) + pw = pap[ipw].w1[:t][:tauzz] + @parallel muladd_tauii!(pw, pv, pac.dmod[:dtM]) + end + end + end +end + +@inbounds @fastmath function add_stress_source!(it::Int64, issp::Int64, iss::Int64, pac::T, pap, activepw, src_flags, sfields=Fields(pac.attrib_mod)) where {T<:P_common{<:FdtdElastic,3}} + if (1 ∈ activepw) + ipw = 1 + sfields = intersect(names(pac.srcwav[ipw][iss].d)[1], sfields) + for sfield in sfields + if (src_flags[ipw]) # add only if src_flags is non-zero + pv = pap[ipw].tauii_buffer # pick stress buffer + pv_vec = view(pv, :) + w = pap[ipw].ss[issp].wavelets[it][sfield] + ssprayw = pap[ipw].ss[issp].ssprayw[sfield] + mul!(pv_vec, ssprayw, w) + pw = pap[ipw].w1[:t][:tauxx] + @parallel muladd_tauii!(pw, pv, pac.dmod[:dtM]) + pw = pap[ipw].w1[:t][:tauyy] + @parallel muladd_tauii!(pw, pv, pac.dmod[:dtM]) + pw = pap[ipw].w1[:t][:tauzz] + @parallel muladd_tauii!(pw, pv, pac.dmod[:dtM]) + end + end + end +end + + + +# # This routine ABSOLUTELY should not allocate any memory, called inside time loop. +@inbounds @fastmath function add_velocity_source!(it::Int64, issp::Int64, iss::Int64, pac, pap, activepw, src_flags, sfields=Fields(pac.attrib_mod)) # adding source to respective sfield at [it] if (1 ∈ activepw) ipw = 1 @@ -72,7 +134,7 @@ end w = pap[ipw].ss[issp].wavelets[it][sfield] ssprayw = pap[ipw].ss[issp].ssprayw[sfield] mul!(pv_vec, ssprayw, w) - @parallel eval(Symbol("muladd_with_density_", sfield))(pw, pv, pac.mod[:rho], pac.fc[:dt]) + @parallel eval(Symbol("muladd_with_density_", sfield, "!"))(pw, pv, pac.mod[:rho], pac.fc[:dt]) end end end @@ -88,22 +150,28 @@ end w = pap[ipw].ss[issp].wavelets[it][sfield] adj_ssprayw = pap[1].ss[issp].rinterpolatew[sfield] mul!(pv_vec, adj_ssprayw, w) - @parallel eval(Symbol("muladd_with_density_", sfield))(pw, pv, pac.mod[:rho], pac.fc[:dt]) + @parallel eval(Symbol("muladd_with_density_", sfield, "!"))(pw, pv, pac.mod[:rho], pac.fc[:dt]) end end end end -# pw = pw + pv * rho -@parallel function muladd_with_density_vx(pw, pv, rho, dt) + +@parallel function muladd_tauii!(pw, pv, dtK) + @all(pw) = @all(pw) + (@all(pv) * @all(dtK)) + return +end + +# # pw = pw + pv * rho +@parallel function muladd_with_density_vx!(pw, pv, rho, dt) @inn(pw) = @inn(pw) + (@inn(pv) / @av_xi(rho) * dt) return end -@parallel function muladd_with_density_vz(pw, pv, rho, dt) +@parallel function muladd_with_density_vz!(pw, pv, rho, dt) @inn(pw) = @inn(pw) + (@inn(pv) / @av_zi(rho) * dt) return end -@parallel function muladd_with_density_vy(pw, pv, rho, dt) +@parallel function muladd_with_density_vy!(pw, pv, rho, dt) @inn(pw) = @inn(pw) + (@inn(pv) / @av_yi(rho) * dt) return end diff --git a/src/fdtd/types.jl b/src/fdtd/types.jl index a728dcd..f54be70 100644 --- a/src/fdtd/types.jl +++ b/src/fdtd/types.jl @@ -83,7 +83,8 @@ mutable struct P_x_worker_x_pw{N,Q<:Data.Array{N}} w1::NamedStack{NamedStack{Q}} # fields stored on GPU or CPU, for both t (current) and tp (previous) time steps w2::NamedStack{NamedStack{Q}} # required for attenuation, where third dimension is nsls (only used for 2D simulation right now) memory_pml::NamedStack{Q} # memory variables for CPML, only derivative fields stored - velocity_buffer::NamedStack{Q} # used for born modeling # only used for 2-D simulation + velocity_buffer::NamedStack{Q} # + tauii_buffer::Q # end # P_x_worker=Vector{P_x_worker_x_pw} @@ -100,6 +101,7 @@ end function reset_w2!(pap::Vector{P_x_worker_x_pw{N,B}}) where {N,B} for pap_x_pw in pap fill!.(pap_x_pw.velocity_buffer, zero(Data.Number)) + fill!(pap_x_pw.tauii_buffer, zero(Data.Number)) map(pap_x_pw.w1) do w fill!.(w, zero(Data.Number)) end diff --git a/src/fwi/fwi.jl b/src/fwi/fwi.jl index 3fdb14d..18f26dc 100644 --- a/src/fwi/fwi.jl +++ b/src/fwi/fwi.jl @@ -49,7 +49,7 @@ function update!(paw::T1, pac::T2) where {T1<:NamedTuple{<:Any,<:Tuple{<:PFdtd,< w = zero(paconv.s) # matrix-free operator A = Conv.operator(paconv, Conv.G()) - IterativeSolvers.lsmr!(w, A, dvec) + IterativeSolvers.lsmr!(w, A, dvec, maxiter=100) copyto!(paconv.s, w) # do conv with estimated w Conv.mod!(paconv, Conv.D()) diff --git a/src/plots.jl b/src/plots.jl index e3fcec1..dd6f20a 100644 --- a/src/plots.jl +++ b/src/plots.jl @@ -29,7 +29,6 @@ @series begin seriestype := :heatmap subplot := j - # seriescolor --> colorschemes[:roma] title --> title1 clims --> Tuple(getfield(medium, field).bounds) seriescolor --> :grays @@ -91,14 +90,14 @@ end layout := (1, length(dat.d)) size --> (length(dat.d) * 300, 300) margin --> 5Measures.mm + dmax = maximum([maximum(dp) for dp in dat.d]) + dmax_clip = dmax - clip_perc * inv(100) * abs(dmax) for (j, field) in enumerate(names(dat.d)[1]) dp = dat[field] - dmax = maximum(abs.(dp)) if (!iszero(dmax)) @series begin dz = dat.grid dx = 1:size(dp, 2) - dmax_clip = dmax - clip_perc * inv(100) * abs(dmax) seriestype := :heatmap framestyle := :grid title --> field