# Seismic modeling and inversion in JUDI.jl

# #Introduction to JUDI

JUDI.jl is a framework for large-scale seismic modeling and inversion and designed to enable rapid translations of algorithms to fast and efficient code that scales to industry-size 3D problems. The focus of the package lies on seismic modeling as well as PDE-constrained optimization such as full-waveform inversion (FWI) and imaging (LS-RTM). Wave equations in JUDI are solved with Devito, a Python domain-specific language for automated finite-difference (FD) computations. JUDI's modeling operators can also be used as layers in (convolutional) neural networks to implement physics-augmented deep learning algorithms. For this, check out JUDI's deep learning extension JUDI4Flux.

The JUDI software is published in Witte, Philipp A., et al. "A large-scale framework for symbolic implementations of seismic inversion algorithms in Julia." Geophysics 84.3 (2019): F57-F71.. For more information on usage of JUDI, you can check the JUDI reference guide.

This tutorial covers the following topics:

- How to set up the geometry/acquisition in a seismic experiment
- How software abstraction in Julia plays in role in the modeling via the abstracted linear operators
- How to set up a mini-cluster to run seismic experiments parallel over shots

`using JUDI, PyPlot, LinearAlgebra`

# #Physical problem setup

## #Grid

JUDI relies on a cartesian grid for modeling and inversion. We start by defining the parameters needed for a cartesian grid:

- A shape
- A grid spacing in each direction
- An origin

```
shape = (201, 201) # Number of gridpoints nx, nz
spacing = (10.0, 10.0) # in meters here
origin = (0.0, 0.0) # in meters as well
```

## #Physical object

JUDI defines a few basic types to handle physical object such as the velocity model. The type `PhyisicalParameter`

is an AbstractVector and behaves as a standard vector.
A `PhysicalParameter`

can be constructed in various ways but always require the origin `o`

and grid spacing `d`

that
cannot be infered from the array.

```
PhysicalParameter(v::Array{vDT}, d, o) where `v` is an n-dimensional array and n=size(v)
PhysicalParameter(n, d, o; vDT=Float32) Creates a zero PhysicalParameter
PhysicalParameter(v::Array{vDT}, A::PhysicalParameter) Creates a PhysicalParameter from the Array `v` with n, d, o from `A`
PhysicalParameter(v::Array{vDT, N}, n::Tuple, d::Tuple, o::Tuple) where `v` is a vector or nd-array that is reshaped into shape `n`
PhysicalParameter(v::vDT, n::Tuple, d::Tuple, o::Tuple) Creates a constant (single number) PhyicalParameter
```

Let's make a simple 3-layer velocity model

```
# Define the velocity (in km/s=m/ms)
vp = 1.5f0 * ones(Float32, shape)
vp[:, 66:end] .= 2.0f0
vp[:, 134:end] .= 2.5f0
# Create a physical parameter
VP = PhysicalParameter(vp, spacing, origin);
```

Let's plot the velocities. Because we adopt a standard cartesian dimension ordering for generality (X, Z) in 2D and (X, Y, Z) in 3D, we plot the transpose of the velocity for proper visualization.

```
figure(figsize=(12, 8));
subplot(121);
imshow(vp', cmap="jet", extent=[0, (shape[1]-1)*spacing[1], (shape[2]-1)*spacing[2], 0]);
xlabel("X [m]");ylabel("Depth [m]");
title("vp");
subplot(122);
imshow(VP', cmap="jet", extent=[0, (shape[1]-1)*spacing[1], (shape[2]-1)*spacing[2], 0]);
xlabel("X [m]");ylabel("Depth [m]");
title("vp as a physical parameter");
```

Because the physical parameter behaves as vector, we can easily perform standard operations on it.

`norm(VP), extrema(VP), 2f0 .* VP, VP .^ 2`

## #Model

JUDI then provide a `Model`

structure that wraps multiple physical parameters together.

`model = Model(shape, spacing, origin, 1f0./vp.^2f0; nb = 80)`

# #Modeling

Now that we have a seismic model, we will generate a few shot records.

## #Acquisition Geometry

The first thing we need is an acquisiton geometry. In JUDI, there are two ways to create a Geometry.

- By hand, as we will show here
- From a SEGY file, as we will show in a follow-up tutorial

We create a split-spread geomtry with sources at the top and receivers at the ocean bottom (top of second layer).

**Note**:

- For 2D simulation (i.e. wave propagation in a 2D plane), JUDI still takes y-coordinate of source/receiver locations but we can just put them to 0 anyway.

```
# Sources position
nsrc = 11
xsrc = range(0f0, (shape[1] -1)*spacing[1], length=nsrc)
ysrc = 0f0 .* xsrc # this a 2D case so we set y to 0
zsrc = 12.5f0*ones(Float32, nsrc);
```

Now this definition creates a single Array of position, which would correspond to a single simultenous source (firing at the same time). Since we are interested in single source experiments here, we convert these position into an Array of Array
of size `nsrc`

where each sub-array is a single source position

`xsrc, ysrc, zsrc = convertToCell.([xsrc, ysrc, zsrc]);`

```
# OBN position
nrec = 101
xrec = range(0f0, (shape[1] -1)*spacing[1], length=nrec)
yrec = 0f0 # this a 2D case so we set y to 0. This can be a single number for receivers
zrec = (66*spacing[1])*ones(Float32, nrec);
```

The last step to be able to create and acquisiton geometry is to define a recording time and sampling rate

```
record_time = 2000f0 # Recording time in ms (since we have m/ms for the velocity)
sampling_rate = 4f0; # Let's use a standard 4ms sampling rate
```

Now we can create the source and receivers geometry

```
src_geom = Geometry(xsrc, ysrc, zsrc; dt=sampling_rate, t=record_time)
# For the receiver geometry, we specify the number of source to tell JUDI to use the same receiver position for all sources
rec_geom = Geometry(xrec, yrec, zrec, dt=sampling_rate, t=record_time, nsrc=nsrc);
```

Let's visualize the geometry onto the model

```
figure();
imshow(vp', cmap="jet", extent=[0, (shape[1]-1)*spacing[1], (shape[2]-1)*spacing[2], 0]);
scatter(xsrc, zsrc, label=:sources);
scatter(xrec, zrec, label="OBN");
xlabel("X [m]");ylabel("Depth [m]");
legend();
title("acquisition geometry")
```

### #Source wavelet

For the source wavelet, we will use a standard Ricker wavelet at 10Hz for this tutorial. In practice, this wavelet would be read from a file or estimated during inversion.

```
f0 = 0.010 # Since we use ms, the frequency is in KHz
wavelet = ricker_wavelet(record_time, sampling_rate, f0);
plot(wavelet)
```

## #judiVector

In order to represent seismic data, JUDI provide the `judiVector`

type. This type wraps a geometry with the seismic data corresponding to it. Note that `judiVector`

works for both shot records at the receiver locations and for the wavelet at the source locations. Let's create one for the source

`q = judiVector(src_geom, wavelet)`

# #Linear operator

Following the math, the seismic data recorded at the receiver locations can be calculated by

Here, $\mathbf{q}$ is the source wavelet set up at the firing locations (defined above). $\mathbf{P}_s^{\top}$ injects $\mathbf{q}$ to the entire space. $\mathbf{A}(\mathbf{m})$ is the discretized wave equation parameterized by squared slowness $\mathbf{m}$. By applying the inverse of $\mathbf{A}(\mathbf{m})$ on source $\mathbf{q}$, we acquire the entire wavefield in the space. Finally, $\mathbf{P}_r$ restricts the wavefield at the receiver locations at the recording time.

Now let's these linear operators, which behave as "matrices".

`Pr = judiProjection(rec_geom) # receiver interpolation`

`Ps = judiProjection(src_geom) # Source interpolation`

`Ainv = judiModeling(model) # Inverse of the disrete ewave equation`

**WARNING**
While these three operator are well defined in JUDI, `judiProjection`

is a no-op operator and cannot be used by itself but only in combination with a `judiModeling`

operator

# #Seismic data generation

With the set-up operators above, we can finally generate synthetic data with simple mat-vec product thanks to the abstraction.

`d_obs = Pr * Ainv * Ps' * q`

Similarly, we can define a modeling operator `F`

that combines the modeling operator with source/receiver restriction operators, as

`F = Pr * Ainv * Ps'`

We can run the same code to generate the data

`d_obs = F * q`

There are 11 sources in this experiment, and the data is generated for each shot independently. If we only want to acquire the data for the 6th source for example, then we can run

`d_obs6 = F[6] * q[6]`

```
data_extent = [xrec[1], xrec[end], 1f-3*record_time, 0]
figure(figsize=(20, 5))
for i=1:5
subplot(1, 5, i)
imshow(d_obs.data[2*i], vmin=-1, vmax=1, cmap="Greys", extent=data_extent, aspect="auto")
xlabel("Receiver position (m)")
ylabel("Recording time (s)")
title("xsrc=$(1f-3xsrc[2*i][1])km")
end
tight_layout()
```

## #Parallelization -- workers and threads

JUDI (based on Devito) uses shared memory OpenMP parallelism for solving PDEs. Here, we show how to build up a small 2-worker cluster by calling `addprocs(2)`

. Because we set the environment variable `ENV["OMP_DISPLAY_ENV"] = "true"`

, we will see the OMP environment printed out on each worker.

We set 4 environment variables related to OpenMP:

`OMP_DISPLAY_ENV`

prints out the OpenMP environment on each worker`OMP_PROC_BIND`

specifies that threads should be bound to physical cores`OMP_NUM_THREADS`

specifies the number of threads per workers is `1/nw`

the number of physical cores`GOMP_CPU_AFFINITY`

specifies which physical cores the threads run on for each worker

If you run the shell command top during execution, you will see 3 julia processes: the main process and 2 workers. The two workers should generally have about 50% of the system, and load average should tend towards the physical number of cores. When running with multiple workers, we gain the parallelism over sources.

```
using Distributed
# Sytem informations
nthread = Sys.CPU_THREADS
nw = 2
ENV["OMP_DISPLAY_ENV"] = "true"
ENV["OMP_PROC_BIND"] = "close"
ENV["OMP_NUM_THREADS"] = "$(div(nthread, nw))"
addprocs(nw; lazy=false)
@show workers()
for k in 1:nworkers()
place1 = (k - 1) * div(nthread,nworkers())
place2 = (k + 0) * div(nthread,nworkers()) - 1
@show place1, place2, div(nthread, nw)
@spawnat workers()[k] ENV["GOMP_CPU_AFFINITY"] = "$(place1)-$(place2)";
end
@everywhere using Distributed
@everywhere using JUDI
```

`d_obs = F * q`

## #Reverse time migration

Reverse time migration (RTM) is a seismic migration method to move the dipping events in the data domain to their supposedly true subsurface positions in the image domain. This is also a linear operation, as shown in our abstraction below.

```
J = judiJacobian(F, q) # forward is linearized born modeling, adjoint is reverse time migration
rtm = J' * d_obs
```

```
figure();
imshow(rtm', cmap="Greys", extent=[0, (shape[1]-1)*spacing[1], (shape[2]-1)*spacing[2], 0]);
xlabel("X [m]");ylabel("Depth [m]");
title("RTM");
```

- Witte, P. A., Louboutin, M., Kukreja, N., Luporini, F., Lange, M., Gorman, G. J., & Herrmann, F. J. (2019). A large-scale framework for symbolic implementations of seismic inversion algorithms in Julia.
*GEOPHYSICS*,*84*(3), F57–F71. 10.1190/geo2018-0174.1