Heterogeneous Map-Reduce: Seismic Imaging Application (Part 2/3)

This article is the second part of the series about the heterogeneous map-reduce approach:
Part 1 – Heterogeneous Map-Reduce: meet the Task System (Part 1/3)
Part 2 – [You are here] – Heterogeneous Map-Reduce: Seismic Imaging Application (Part 2/3)
Part 3 – Heterogeneous Map-Reduce: Scientific Visualisation (Part 3/3)

In the previous part we have described the heterogeneous map reduce framework. Here, we will start with an example from a seismic imaging application.

Problem setting

The oil and gas industry makes use of computational intensive algorithms such as reverse-time migration and full waveform inversion to provide an image of the subsurface. The image is obtained by sending a wave energy into the subsurface and recording the signal required for a seismic wave to reflect back to the surface from the interfaces with different physical properties. A seismic wave is usually generated by shots at known frequencies, placed close to the surface on land or to the water surface in the sea. Returning waves are usually recorded in time by hydrophones in the marine environment or by geophones during land acquisition. The goal of seismic imaging is to transform the seismograms to a spatial image of the subsurface.

Migration algorithms produce an image of the subsurface given seismic data measured at the surface. In particular, pre-stack depth migration produces the depth locations of reflectors by mapping seismic data from the time domain to the depth domain, assuming that a sufficiently accurate velocity model is available. The classic imaging principle is based on the correlation of the forward propagated wavefield from a source and a backward propagated wavefield from the receivers. In frequency domain the image is calculated as follows:

\text{Image} (x,y,z) = \sum_{\text{shots}}\sum_{\omega} W_{\text{source}}^*(x,y,z,\omega) W_{receivers}(x,y,z,\omega)

where W_{\text{source}}(x,y,z,\omega)  denotes the wavefield propagated from the source and W_{\text{receivers}}(x,y,z,\omega)  from the receivers respectively, \omega denotes the frequency. That means for every shot and every frequency we need to simulate the wave propagation twice.

Levels of parallelism

Level 1: The highest level of parallelization for frequency-domain migration is over the shots. Each shot is treated independently. We assume that the migration volume for one shot is computed on one compute node connected to none, one or more GPUs.

Level 2: The next level of parallelism involves the frequencies. For each frequency, a linear system of equations needs to be solved.

Level 3: The third level of parallelism includes matrix decomposition, where the matrix for the linear system of equations is decomposed into subsets of rows that fit on a single GPU.

Level 4: The last level of parallelism for migration in frequency domain is parallelization of matrix-vector multiplications (MVMs) and vector-vector operations.

Task system (Map-Reduce)

Here we describe how the task system works for the migration in frequency domain.

Read more about task system in the previous post: Heterogeneous Map-Reduce: Meet the Task System

 

Mostly we will use the first level of the parallelism. The server or ‘master node’ creates one task per source. Each task is added to a “Available” queue. When a client requests a task, a given task is moved from the queue to the “Running” list.

As we have seen earlier, the migration algorithm consists of forward and backward modelling in frequency domain for each source. Therefore, the second level of parallelism for migration consist of parallelization over all frequencies for each source.

Let’s assume that we have N_s sources and N_\omega frequencies. Then, one task consists of computations of one frequency \omega_{s_i} for a given source {s_i}, i=1,...,N_s. In total, we have N_\omega\cdot N_s tasks.

For each frequency, a linear system of equations needs to be solved.  The matrix size and memory requirements are the same for each frequency, but the lower frequencies require less compute time than the higher ones. Here, we assume that one frequency for one source in the frequency domain fits in one compute node. At this point, the auto load-balancing of the tasks comes into play. There is no need to know beforehand how to distribute the shots over the compute nodes.

If a compute node is connected to one or more GPUs, we can make use of the third level of parallelism and decompose the matrix across GPU(s). However, this is done within a task.

Migration movie

The movie starts with the velocity model. We run migration on the SEG/EAGE Overtrust model, which represents an acoustic constant-density medium with complex, layered structures and faults. Then the camera rotates and it shows the migrated image, which appears to be empty. The migrated data from each shots is added and shown in the movie as it is being received by the server (‘Reduce’ part of the algorithm).  We can see how layers become visible. At the end we see the final image that we obtained using migration in frequency domain.

More details on the simulation

The volume has a size of 1000×1000×620 m3. The problem is discretized on a grid with 301×301×187 points and a spacing of 3.33 m in each coordinate direction.

The discretization for migration in frequency domain is 2nd-order in space. A Ricker wavelet with a peak frequency of 15 Hz is chosen for the source and the maximum frequency in this experiment is 30 Hz. Note that by reducing the maximum frequency, we can increase the grid spacing. For instance, by choosing a maximum frequency of 8 Hz, the grid spacing can be chosen as 25 m in each direction. The line of sources is located at a depth of 10 m and is equally spaced with an interval of 18.367 m 4 in the x-direction.

The receivers are equally distributed in the two horizontal directions with the same spacing as the sources, at the depth of 20 m. The sampling interval for the modelled seismic data is 4 ms. The maximum simulation time is 0.5 s. For migration in the frequency domain we chose a frequency interval of 2 Hz.

What is next?

In this post we showed how to use the heterogeneous map-reduce framework and the task system on an application from seismic imaging : migration in frequency domain.

In the next post of this series we are going to have a look at scientific visualisation.


Get EZNumeric’s future articles in your inbox:

2 comments on “Heterogeneous Map-Reduce: Seismic Imaging Application (Part 2/3)Add yours →

Leave a Reply

Your email address will not be published. Required fields are marked *