# SeisNoise.jl GPU Computing Tutorial

* Abstract *We introduce SeisNoise.jl, a library for high-performance ambient seismic noise cross-correlation, written entirely in the computing language Julia. Julia is a new language, with syntax and a learning curve similar to MATLAB, R, or Python and performance close to Fortran or C. SeisNoise.jl is compatible with high-performance computing resources using both the Central Processing Unit (CPU) and the Graphics Processing Unit (GPU). SeisNoise.jl is a modular toolbox, giving researchers common tools and data structures to design custom ambient seismic cross-correlation workflows in Julia.

### Check out the Paper in SRL!

*SeisNoise.jl: Ambient Seismic Noise Cross-Correlation on the CPU and GPU in Julia*

Tim Clementsinline_formula not implemented, and Marine A. Denolleinline_formula not implemented

inline_formula not implementedDepartment of Earth and Planetary Sciences, Harvard University, Cambridge, MA

This is a tutorial to accompany paper in *Seismological Research Letters*

### Why Julia?

### Julia is fast!

Julia was designed from the beginning for high performance. Julia programs compile to efficient native code for multiple platforms via LLVM.

### Dynamic

Julia is dynamically typed, feels like a scripting language, and has good support for interactive use.

### Optionally typed

Julia has a rich language of descriptive datatypes, and type declarations can be used to clarify and solidify programs.

### General

Julia uses multiple dispatch as a paradigm, making it easy to express many object-oriented and functional programming patterns. It provides asynchronous I/O, debugging, logging, profiling, a package manager, and more.

### Easy to use

Julia has high-level syntax, making it an accessible language for programmers from any background or experience level. Browse the Julia microbenchmarks to get a feel for the language.

### Open source

Julia is provided under the MIT license, free for everyone to use. All source code is publicly viewable on GitHub.

### Learning Julia

We recommend the learning Julia page for resources on learning Julia.

## Installation

SeisNoise.jl is a registered Julia package and can be installed using the Julia package manager. Here we'll install SeisNoise.jl and a couple other packages, such as SeisIO.jl for loading/downloading data and CUDA.jl for using the GPU.

`using Pkg `

`Pkg.add("SeisIO")`

`Pkg.add("CUDA")`

`Pkg.add("BenchmarkTools")`

`Pkg.add(PackageSpec(url="https://github.com/tclements/SeisNoise.jl", rev="master")) # add the most recent version from github`

Before we can use a package in Julia, we need to import it with the `using`

syntax. We'll load the following modules

`Dates`

: For time manipulation`CUDA`

: Interface to GPU computations`SeisIO`

: For loading/downloading/pre-processing seismic data`SeisNoise`

: Ambient Noise Cross-Correlation

The first time these modules are loaded will take some time, as they need to be pre-compiled...

`using BenchmarkTools, CUDA, Dates, SeisIO, SeisNoise `

## Brief Primer on Ambient Noise Cross-Correlation

Given two seismometers, inline_formula not implemented and inline_formula not implemented, placed at some distance apart on the surface, the stations will record ground motion as a function of time: inline_formula not implemented and inline_formula not implemented. Over long periods of time, the cross-correlation of ground motions of inline_formula not implemented and inline_formula not implemented

inline_formula not implemented

yields a band-limited approximation of the elastodynamic Green's function, or impulse response, between inline_formula not implemented and inline_formula not implemented. We call inline_formula not implemented the time-domain cross-correlation function. Using the convolution theorem, we can re-write cross-correlation in the frequency domain:

inline_formula not implemented

where inline_formula not implemented is the frequency-domain cross-correlation, inline_formula not implemented is the Fourier transform of inline_formula not implemented and inline_formula not implemented is the Fourier transform of inline_formula not implemented, and inline_formula not implemented denotes the complex conjugate. We can go retrieve the time-domain cross-correlation from the frequency-domain cross-correlation using the inverse Fourier transform:

inline_formula not implemented

SeisNoise.jl computes cross-correlations using inline_formula not implemented and inline_formula not implemented rather than inline_formula not implemented due to the efficiency of the Fast Fourier Transform (FFT). We'll show in the following sections how to implement this with SeisNoise.jl.

## SeisNoise Structure

Processing in SeisNoise.jl follows a flow of data through custom structures for ambient noise cross-correlation.

We use the SeisIO.jl module to load seismic data. SeisIO.jl can read data in the following formats:

ASDF, GeoCSV, QuakeML, SAC, SEED, SEG Y, SLIST, SUDS, StationXML, Win32

SeisIO.jl provides two types of structures for handling data:

`SeisChannel`

: stores single-channel geophysical data

`SeisData`

: stores multi-channel geophysical

SeisNoise.jl provides three custom structures for computing ambient noise cross-correlations:

`RawData`

: stores ambient seismic noise data in short, overlapping time windows -> inline_formula not implemented and inline_formula not implemented, as above`FFTData`

: stores Fourier transforms of these time windows -> inline_formula not implemented and inline_formula not implemented, as above`CorrData`

: stores the corresponding NCFs -> inline_formula not implemented, as above

### Cross-Correlation Pre-Processing

SeisNoise.jl uses SeisIO.jl for reading data and pre-processing. Here we'll use the Southern California Earthquake Data Center (SCEDC) Open Dataset on Amazon Web Services (AWS) to test some cross-correlation with SeisNoise.jl.

First let's load some data using the SeisIO `get_data`

function:

`S = read_data("/scedc-pds/continuous_waveforms/2019/2019_031/CISDD__BHZ___2019031.ms")`

S holds metadata, processing notes, timing information and data for each channel. Below is a simple SeisIO.jl workflow to downsample `S`

from 40 to 20 Hz:

`ungap!(S) # remove gaps`

`detrend!(S) # remove mean and trend`

`taper!(S) # taper ends`

`resample!(S,fs=20.) # downsample from 40 to 20 [Hz]`

Functions in Julia that end in ! (e.g. detrend!(S) as above), by convention modify their arguments in-place, while functions without !’s (e.g. detrend(S)) allocate new arrays or structures.

`resample(S,fs=10.) # allocates new SeisData structure `

`S`

Note that the call to `resample(S,fs=10.)`

does not change `S`

's sampling rate to 10 Hz, but rather allocates a new `SeisData`

, leaving `S`

unchanged.

### Pre-Processing with RawData

SeisNoise.jl provides the `RawData `

structure for pre-processing windows of ambient noise. The advantage of `RawData`

over `SeisData`

is data can be processed as a matrix. This is especially helpful when processing on the GPU. To create a `RawData`

structure, we need to specify:

`cc_len`

: length of each time window in seconds`cc_step`

: step between time windows in seconds

RawData structures can be created from either SeisData or SeisChannel:

`cc_len = 1800 # time window length [s] = 30 minutes `

`cc_step = 450 # step between time windows [s] = 7.5 minutes or 75% overlap `

`R = RawData(S,cc_len,cc_step) # RawData struct`

Ambient noise data is accessible via the .x field, while .t gives the start time of each window in UNIX time. Here, we'll filter the RawData between 0.1 and 0.3 Hz, apply clipping and then whiten the data between 0.12 and 0.28 Hz.

`freqmin,freqmax = 0.1,0.3 # min/max freqs [Hz]`

`detrend!(R) # detrend each window `

`taper!(R) # taper `

`bandpass!(R,freqmin,freqmax)`

`clip!(R,3) # clip data at 3x RMS of each window`

`whitemin, whitemax = 0.12, 0.28 # whiten min/max [Hz]`

`whiten!(R,whitemin,whitemax) # apply whitening`

### Fourier Transforms with FFTData

Once the `RawData `

is pre-processed, we take the Fast Fourier Transform (fft) to convert to the frequency domain. In SeisNoise.jl the `rfft`

function computes the real Fourier transform of the data in a `RawData`

structure and returns a `FFTData`

structure. We use the real fft because it is 2-3x faster than a regular fft ambient noise data.

`FFT = rfft(R)`

`FFTData `

allow for spectral smoothing or whitening before cross-correlating

`whiten!(FFT,whitemin,whitemax)`

### Cross-correlating with SeisNoise

Cross-correlating with Seisnoise is done with the `correlate`

function. `correlate `

accepts two `FFTData`

structures and a `maxlag`

value - the maximum number of seconds of coda to return, as input. Here we'll input the same FFTData twice, giving an autocorrelation.

`maxlag = 60. # maximum lag time in correlation`

`C = correlate(FFT,FFT,maxlag)`

Note that correlation data is stored in the .corr field. CorrData have a few new parameters, namely

comp: cross-correlation component (e.g. "ZZ","EN","RT")

rotated: true/false for if the CorrData has been rotated to ZRT reference frame

corr_type: Type of correlation: `cross-correlation`, `coherence` or

`deconv`.

dist: Distance between stations in km

azi: Azimuth from first station to second station

baz: Backazimuth from first station to second station

maxlag: The maximum number of seconds of coda

### Post-Processing

The most common post-processing is stacking. The SeisNoise.jl `stack`

function can stack arbitrary time periods. The default is to stack by day:

`stack(C)`

But you can specify the stack interval from seconds, hours, days, months, years, etc..

`stack(C,interval=Minute(30))`

`stack(C,interval=Hour(2))`

## SeisNoise.jl with GPUs

Why Graphical Processing Units (GPUs)?

GPU have many cores vs few cores on a PC. GPU throughput is often an order of magnitude greater than CPU throughput. This is good for linear algebra + computationally intensive tasks aka Ambient Noise Cross Correlation.

### GPU Computing in Julia

We use the CUDA.jl package. CUDA.jl makes it easy to move data onto NVIDIA GPUs using the `cu`

function

`A = cu(rand(Float32,10,10))`

Array `A `

is on the GPU. We can do most array operations with `A`

, such as

`A + A # add A to itself `

`A .+ 2 # add 2 to each element`

`A .- 4 # subtract 4 from each element`

`A * A # multiply A with itself`

`A .* A # square each element`

Let's look at the time difference between matrix multiplication on the CPU and GPU

`B = rand(Float32,1000,1000) # this array is on the CPU `

`C = cu(B); # on the GPU `

`b = B * B # benchmark on the CPU`

`c = CUDA. C * C # benchmark on the GPU `

`ratio(minimum(b),minimum(c))`

Matrix multiplication is ~14x faster on the GPU! Let's leverage this power for ambient noise cross-correlation in the next section.

### GPUs + SeisNoise

GPU computing with SeisNoise.jl is relatively easy - there's no need to write any GPU code. Use the `gpu`

function to move data stored in RAM on the CPU to the GPU

`RGPU = R |> gpu `

The RawData structure `RGPU`

is now on the GPU! (in actuality, only R.x) is on the GPU. We can apply the same processing as before but on the GPU

`detrend!(RGPU) # detrend each window `

`taper!(RGPU) # taper `

`bandpass!(RGPU,freqmin,freqmax)`

`clip!(RGPU,3) # clip data at 3x RMS of each window`

`whitemin, whitemax = 0.12, 0.28 # whiten min/max [Hz]`

`whiten!(RGPU,whitemin,whitemax) # apply whitening`

Let's look at the performance of detrending on the CPU and GPU

`dCPU = detrend!(R)`

`dGPU = detrend!(RGPU)`

`ratio(minimum(dCPU),minimum(dGPU))`

It's 23x faster! Taking an fft of a `RawData`

works similarly on the GPU

`FFTGPU = rfft(RGPU)`

Once we have ffts on the GPU, we can cross-correlate on the GPU.

`CGPU = correlate(FFTGPU,FFTGPU,maxlag)`

## Full GPU Cross-correlation workflow

Here is a complete cross-correlation routine on the GPU. Most of the time is spent loading data from Amazon Web Services. The only difference to a CPU-based processing routine is calling `|> gpu`

to send waveform data to the GPU and `|> cpu`

to send correlation back to the CPU.

`# load SCEDC data from AWS `

`S1 = S = read_data("/scedc-pds/continuous_waveforms/2019/2019_031/CISDD__BHZ___2019031.ms")`

`S2 = S = read_data("/scedc-pds/continuous_waveforms/2019/2019_031/CIPER__BHZ___2019031.ms")`

`# convert SeisData to RawData and send to GPU`

`R1 = RawData(S1, cc_len, cc_step) |> gpu`

`R2 = RawData(S2, cc_len, cc_step) |> gpu`

`R = [R1,R2]`

`# preprocess on the GPU`

`detrend!.(R)`

`taper!.(R)`

`bandpass!.(R,freqmin,freqmax,zerophase=true)`

`# Real FFT on GPU`

`FFT = rfft.(R)`

`whiten!.(FFT,freqmin,freqmax)`

`# compute correlation and send to cpu`

`C = correlate(FFT[1],FFT[2],maxlag) |> cpu`

Now let's plot the correlation:

`corrplot(C)`

### Extending SeisNoise.jl

Let's say you want to make your own function to work with RawData, FFTData, or CorrData structures. We recommend using multiple dispatch to write kernel functions for computations on Arrays and then writing higher-level functions that take RawData, FFTData, or CorrData as input, then call kernel functions. Here's an example to double the amplitude of data:

`# here is a kernel function `

`function double!(A::AbstractArray{T}) where T <: Real# here T is the element type `

` A .*= T(2)`

` return nothing `

`end`

`# here is function for RawData`

`double!(R::RawData) = double!(R.x)`

`# here is function for FFTData`

`double!(F::FFTData) = double!(F.fft)`

`# and here is CorrData`

`double!(C::CorrData) = double!(C.corr)`

There are now 4 different versions of `double!`

. Julia will pick the correct one depending on input:

`A = collect(1:4)`

`double!(A)`

`A`

And now we'll try it on a RawData structure

`R = RawData() # empty structure`

`R.x = reshape(collect(1.:9),3,3) |> Array`

`double!(R)`

`R.x`

If you're interested in what's going on here, we can use the `@code_lowered`

to show the byte code involved

`double!(R) `

`double!(R)`

first gets `R.x`

, then calls the kernel version of `double!`

, as shown below

`double!(A) `

which is just a broadcast operation. If you've made it this far and are interested in contributing, take a look at the SeisNoise.jl source code on Github. Happy correlating!