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 manipulationCUDA
: Interface to GPU computationsSeisIO
: For loading/downloading/pre-processing seismic dataSeisNoise
: 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 aboveFFTData
: stores Fourier transforms of these time windows -> inline_formula not implemented and inline_formula not implemented, as aboveCorrData
: 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 secondscc_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!