Skip to content

Commit

Permalink
SeisForwExpt docs updated
Browse files Browse the repository at this point in the history
  • Loading branch information
pawbz committed Dec 16, 2019
1 parent 95d26e1 commit e1eb677
Show file tree
Hide file tree
Showing 19 changed files with 391 additions and 406 deletions.
6 changes: 3 additions & 3 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,9 @@ makedocs(
"SrcWav" => "generated/srcwav/doc.md",
"Data" => "generated/data/doc.md",
"SeisForwExpt" => Any[
"Introduction" => "generated/fdtd/doc.md",
"Basic usage" => "generated/fdtd/reuse_expt.md",
"Generate snaps" => "generated/fdtd/create_snaps.md",
"Description" => "generated/fdtd/doc.md",
"Tutorial I" => "generated/fdtd/create_snaps.md",
"Tutorial II" => "generated/fdtd/reuse_expt.md",
],
"PoissonExpt" => Any[
"Introduction" => "generated/Poisson/doc.md",
Expand Down
87 changes: 50 additions & 37 deletions docs/src/generated/fdtd/create_snaps.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,70 +2,83 @@
EditURL = "<unknown>/fdtd/create_snaps.jl"
```

We aim to save the snapshots, at given time steps in `SeisForwExpt`.

### Loading some packages

```@example create_snaps
using GeoPhyInv
using Statistics
using Plots; gr();
nothing #hide
```

### Setting up the variables necessary to create the `Expt`
### Medium

```@example create_snaps
model=Medium(:acou_homo1); # load a simple homogeneous acoustic model from the gallery
update!(model, [:vp,:rho], randn_perc=0.1); # add some random noise to the model
ageom=AGeom(model.mgrid,:xwell, SSrcs(2)); # load a simple acquisition ageometry using `mgrid` of the seismic model
tgrid = range(0.0,stop=2.0,length=2000) # generate a time grid
wav = ricker(10.0, tgrid, tpeak=0.25,); # ricker wavelet
medium=Medium(:acou_homo1); # load a simple homogeneous acoustic medium from the gallery
update!(medium, [:vp,:rho], randn_perc=5); # add some random noise to the medium
nothing #hide
```

srcwav = SrcWav(tgrid, ageom, [:P])
update!(srcwav, [:P], wav)
@info "We are ready for the modeling."
### AGeom

```@example create_snaps
ageom=AGeom(medium.mgrid,:xwell, SSrcs(2)); # load a simple acquisition using `mgrid` of the medium
nothing #hide
```

### Final step
### SrcWav

One can plot the model, source and receivers using these commands:
`using Plots;`
`p1=JP.seismic(model);`
`JP.ageom!(ageom);`
`plot(p1);`
```@example create_snaps
tgrid = range(0.0,stop=2.0,length=2000); # generate a time grid
wav = ricker(10.0, tgrid, tpeak=0.25,); # ricker wavelet
srcwav = SrcWav(tgrid, ageom, [:P]);
update!(srcwav, [:P], wav);
nothing #hide
```

### Plotting #1

```@example create_snaps
p1=heatmap(medium, :vp)
scatter!(ageom, SSrcs())
scatter!(ageom, Recs())
plot(p1)
```

### SeisForwExpt
Now we have all the required variables to create `SeisForwExpt` object and
prepare the forward modelling.
While creating, we switched the `snaps_flag` on, and instructed recording field at
prepare the forward modeling.
While creating, we switched the `snaps_flag` on, and instruct recording field at
`tsnaps`.
Once the `Expt` object is created, do the modelling "without approximately any
memory allocations" using `mod!`
Once the `Expt` object is created, do the modeling "without approximately any
memory allocations" using `update!`

```@example create_snaps
paE=SeisForwExpt(Fdtd(),model=model,
ageom=[ageom], srcwav=[srcwav],
snaps_flag=true,
tsnaps=[0.3, 0.4, 0.5],
tgridmod=tgrid, verbose=true);
pa=SeisForwExpt(Fdtd(),medium=medium, ageom=ageom, srcwav=srcwav,
snaps_flag=true, tsnaps=[0.3, 0.4, 0.5], tgrid=tgrid, verbose=true);
nothing #hide
```

### Modeling

@time update!(paE);
```@example create_snaps
@time update!(pa);
nothing #hide
```

### Extracting snaps from `Expt`
### Extracting snaps

```@example create_snaps
snaps=paE[:snaps,1]; # extracting snaps of the first supersource
snaps=paE[:snaps,2]; # second supersource
@info string("The dimensions of the snaps are (nz,nx,nt)=", size(snaps))
snaps=pa[:snaps,1]; # extracting snaps of the first supersource
snaps=pa[:snaps,2]; # second supersource
nothing #hide
```

We can now plot snapshots using these commands:
### Plotting #2

```@example create_snaps
p=[]
#Gadfly.set_default_plot_size(20cm, 5cm)
for ii in 1:3
# push!(p, spy(snaps[:,:,ii], Guide.xlabel("x"), Guide.ylabel("z")))
push!(p,heatmap(medium.mgrid[2], medium.mgrid[1], snaps[:,:,ii], xlabel="x", ylabel="z"))
end
#hstack(p...)
plot(p..., layout=(1,3), size=(1100,250))
```

22 changes: 17 additions & 5 deletions docs/src/generated/fdtd/doc.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,23 @@
EditURL = "<unknown>/fdtd/doc.jl"
```

A forward experiment, where the seismic data are generated
using some models and acquisition parameters from our gallery.
Forward modeling consists of a finite-difference simulation of the acoustic wave-equation.

```@example doc
b=4
using BenchmarkTools
using GeoPhyInv
using Test
```

# Intro

```@docs
SeisForwExpt
Base.getindex(::GeoPhyInv.PFdtd, ::Symbol, ::Int)
```

# Methods
```@docs
update!(::GeoPhyInv.PFdtd)
update!(::GeoPhyInv.PFdtd, ::Medium)
update!(::GeoPhyInv.PFdtd, ::SrcWav, ::Any)
```

121 changes: 47 additions & 74 deletions docs/src/generated/fdtd/reuse_expt.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,141 +2,114 @@
EditURL = "<unknown>/fdtd/reuse_expt.jl"
```

This tutorial first creates a variable `pa`, i.e. allocates
necessary memory to perform `SeisForwExpt`.
Then, we perform forward modeling using an in-place
function `mod!`. Finally, we will update `pa` with a different
subsurface model and re-run the modeling task with no additional memory
allocation.
The ability to iteratively run the forward modeling task on
various subsurface models is necessary while implementing inversion
algorithms.

### Load packages

```@example reuse_expt
using GeoPhyInv
using Statistics
using Plots
using Plots; gr();
using Test
```

### Setup a spatial grid
This tutorial demonstrates how an instance of `SeisForwExpt` can be used iteratively
to perform forward modeling using `update!` for
various bundles of medium parameters.

### Grid

```@example reuse_expt
zgrid=range(-1000.0,stop=1000.0,length=201)
xgrid=range(-1000.0,stop=1000.0,length=201)
mgrid = [zgrid, xgrid]
@info string("spatial sampling intervals (dz,dx)=", step.(mgrid))
mgrid = [zgrid, xgrid];
nothing #hide
```

### Allocate a `Seismic` model, and adjust bounds
### Medium #1

```@example reuse_expt
vpb = [1500., 3500.] # bounds for vp
rhob = [1500., 3500.] # density bounds
model=Medium(mgrid)
update!(model, [:vp,:rho], [vpb,rhob])
fill!(model)
@info "a seismic model is created"
medium=Medium(mgrid);
update!(medium, [:vp,:rho], [vpb,rhob]);
fill!(medium);
update!(medium, [:vp,:rho], randn_perc=5); # add some random noise
nothing #hide
```

### Add some noise to the model (optional)
### Medium #2

```@example reuse_expt
update!(model, [:vp,:rho], randn_perc=0.01); # add some random noise
nothing #hide
medium_new=similar(medium)
update!(medium_new, [:vp,:rho], randn_perc=5.) # add some noise
update!(medium_new, [:vp], constant_pert=500.) # perturb velocity
```

### A surface acquisition ageometry
### AGeom

```@example reuse_expt
ageom=AGeom(model.mgrid,:surf, SSrcs(3), Recs(30));
ageom=AGeom(medium.mgrid,:surf, SSrcs(3), Recs(30)); # surface seismic
nothing #hide
```

### Plot the model and source, receivers
### Plotting #1

```@example reuse_expt
p1=heatmap(model, :vp)
p1=heatmap(medium, :vp)
scatter!(ageom, SSrcs())
scatter!(ageom, Recs())
plot(p1)
p2=heatmap(medium_new, :vp)
scatter!(ageom, SSrcs())
scatter!(ageom, Recs())
plot(p1, p2, size=(800,300))
```

### Generate a temporal grid
### SrcWav

```@example reuse_expt
tgrid = range(0.0,stop=2.0,length=1000)
tgrid = range(0.0,stop=1.0,length=1000) # generate a temporal grid
wav = ricker(10.0, tgrid, tpeak=0.25,); # Choose a source wavelet
srcwav = SrcWav(tgrid, ageom, [:P]) # initialize
update!(srcwav, [:P], wav) # distribute to all supersources
```

### Choose a source wavelet
### SeisForwExpt

```@example reuse_expt
wav = ricker(10.0, tgrid, tpeak=0.25,);
pa=SeisForwExpt(Fdtd(),medium=medium, ageom=ageom, srcwav=srcwav, tgrid=tgrid, verbose=true);
nothing #hide
```

### Distribute the same source wavelet to all the supsersources

```@example reuse_expt
srcwav = SrcWav(tgrid, ageom, [:P])
update!(srcwav, [:P], wav)
```

create `Fdtd.Param` object to prepare forward modelling
* npw corresponds to the number of independently propagating wavefields (1 in most cases)
Once the `Param` object is created, do the modelling "without any memory allocations" using `mod!`
### Modeling #1

```@example reuse_expt
pa=SeisForwExpt(Fdtd(),npw=1,model=model,
ageom=[ageom], srcwav=[srcwav],
sflags=[2], rflags=[1],
tgridmod=tgrid, verbose=true);
@time update!(pa);
d1=copy(pa[:data][:P])
p1=heatmap(pa[:data]);
nothing #hide
```

plot a record after modelling
### Update medium in `pa` without memory allocation

```@example reuse_expt
pdata=heatmap(pa[:data])
plot(pdata)
update!(pa, medium_new)
```

create new seismic model
### Modeling #2 and plot

```@example reuse_expt
model_new=Medium(:acou_homo1) # prepare another model
update!(model_new, [:vp,:rho], randn_perc=0.01)
update!(model_new, [:vp,:rho], constant_pert=0.03) # perturb the model
p2=heatmap(model_new, :vp) # plot new model
scatter!(ageom, SSrcs())
scatter!(ageom, Recs())
plot(p2)
```

Now, we the change the model in the `Param` object without memory allocation
This routine can be used during FWI,
where medium parameters are itertively updated in the same `Fdtd.Param` object

```@example reuse_expt
update!(pa, model_new)
@time update!(pa);
d2=copy(pa[:data][:P])
p2=heatmap(pa[:data]);
nothing #hide
```

run modelling now and plot data again
Test

```@example reuse_expt
@time update!(pa);
nothing #hide
@test d1≠d2
```

plot a record after modelling
### Plotting

```@example reuse_expt
heatmap!(pdata, pa[:data])
plot(pdata)
plot(p1,p2, size=(500, 300))
```

6 changes: 4 additions & 2 deletions docs/src/index.md
Original file line number Diff line number Diff line change
@@ -1,11 +1,13 @@
# The Expt Datatypes
# The Expt Types

The methods in this package numerically solve some
differential equations commonly faced in geophysical inverse problems.
The functionality of this package revolves around the mutable `Expt` types.
Julia's multiple dispatch is used to overload `Base` methods whenever possible.
Which means, if one is familiar with the `Base` methods of Julia, then (almost) no additional syntax is required
to use this package.
to use this package. Easy!


While performing a given experiment,
firstly, most of the memory necessary
is allocated while creating the `Expt` variables.
Expand Down
10 changes: 5 additions & 5 deletions src/Born/Born.jl
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ end
function mod(;
vp0::Float64=2000.0,
rho0::Float64=2000.0,
model_pert::Medium=nothing,
medium_pert::Medium=nothing,
born_flag::Bool=false,
tgridmod::StepRangeLen=nothing,
tgrid::StepRangeLen = tgridmod,
Expand All @@ -42,15 +42,15 @@ function mod(;


if(born_flag)
mesh_x=model_pert.mgrid[2]
mesh_z=model_pert.mgrid[1]
mesh_x=medium_pert.mgrid[2]
mesh_z=medium_pert.mgrid[1]
nz=length(mesh_z)
nx=length(mesh_x)
δx = step(mesh_x)
δz = step(mesh_z)

δmodtt = model_pert[:KI] - (vp0 * vp0 * rho0)^(-1)
δmodrr = model_pert[:rhoI] - (rho0)^(-1)
δmodtt = medium_pert[:KI] - (vp0 * vp0 * rho0)^(-1)
δmodrr = medium_pert[:rhoI] - (rho0)^(-1)
end


Expand Down
Loading

0 comments on commit e1eb677

Please sign in to comment.