Parallel Simulation of 3D Wave Propagation by Domain Decomposition

doi:10.4236/jamp.2013.14002

Paper Menu >>

Journal Menu >>

Journal of Applied Mathematics and Physics, 2013, 1, 6-11

http://dx.doi.org/10.4236/jamp.2013.14002 Published Online October 2013 (http://www.scirp.org/journal/jamp)

Parall el Simul at ion of 3D Wave Propagati on by Doma in

Decomposition

Galina Reshetova1, Vladimir Tcheverda2, Dmitry Vishnevsky2

1The Institute of Computational Mathematics and Mathematical Geophysics, Siberian Branch of the RAS, Russia

2The Institute of Petroleum Geology and Geophysics, Siberian Branch of the RAS, Russia

Email: [email protected]

Received July 2013

ABSTRACT

In order to perform large scale numerical simulation of wave propagation in 3D heterogeneous multiscale viscoelastic

media, Finite Difference technique and its parallel implementation based on domain decomposition is used. A couple of

typical statements of borehole geophysics are dealt with—sonic log and cross well measurements. Both of them are

essentially multiscales, which claims to take into account heterogeneities of very different sizes in order to provide re-

liable results of simulations . Locally refined spatial grids help us to avoid the use of redundantly tiny grid cells in a tar-

get area, but cause some troubles with uniform load of Processor Units involved in computations. We present r esults of

scalability tests together with results of numerical simulations for both statements performed for some realistic models.

Keywords: Seismic Wave propagation; Sonic Log; Numerical Simulation; Domain Decomposition

1. Introduction

The most effective way to improve resolving ability of

any wave ima ges is to increase dominant frequency of a

sounding pulse. But Earth media attenuate and disperse

propagating waves. Both these effects are often quanti-

fied by the Quality Factor Q. This factor describes rela-

tive dissipation of seismic energy per unit volume per

unit cycle. If attenua tion is not too strong, Quality Factor

can be treated as a number of wavelengths a wave can

propagate through a medium before its amplitude was

decreased in

time s.

Both fields and laboratory experiments prove that this

parameter can be treated as independent on time fre-

quency for rather wide frequency range [5]. Therefore

the higher is the dominant frequency of a source pulse,

the shorter is the distance with reliable level of sig-

nal-to-noise ratio. So, in order to get an image with high

resolution it is necessary to place acquisition system as

close to the target object as possible. The only way to do

this is to place sources and/or receivers within boreholes

drilled in the vicinity of the targ et object. I n its own turn

presence of a well filled with a drilling mud brings es-

sential peculiarities to wave fields and should be taken

into account in model description. In its own turn, this

claims use of locally refined grids in order to catch hete-

rogeneities of the smallest scale only.

The paper deals with numerical simulation of waves

propagation within viscoelastic media for the following

borehole ba s ed geophysical methods:

a) Sonic Log—sources and receivers are within the

same borehole and the main task is monitoring of casing

pipe integrity and recovery of near well-bore vicinity.

Multiscale nature of the problem becomes apparent in the

presence of heterogeneities of at least two different

scales-distance source/receiver and borehole radius. If

one deals with a cased borehole there is third scale-

structure of the casing. Fourth scale can be introduced by

a medium;

b) Cros s-well Tomography—sources and receivers

are placed within adjacent boreholes encircling a target

object. This problem possesses two extremely different

scales-borehole diameter and distance between sources

and receivers.

2. Viscoelastic Media

Any process of wave propagation within linear elastic

medium is governed by t w o groups of eq ua tions:

1) Motion equations (Newton’s law);

2) State equation that connects stress and strain tensors

(Hook’s law).

In reality rocks possess intrinsic attenuation produced

by their “memo ry”: stress state at some instant t is de-

termined by the “history” of strains. We will suppose that

all processes are causal, that is the current state of the

medium does not depend on the future by no means. So,

in the most general form the Hook’s law for these mate-

G. RESHETOVA ET AL.

rials can be written down in the following way:

( )( )

(,) (,0)(,)

ijijkl kl

tkl

ijkl

xtG xxt

Gx d

σε

ετ

ττ

= +

∂−

+∂

∫

For isotropic viscoelastic materials relaxation tensor

simplifies to the follo wing one:

Λ(,) 2(,)

ijklij klikjl

GxtM xt

δδ δδ

= +

Numerical resolution of integro-differential equations

(1) is very troublesome. The most convenient way is to

represent state equation in a differential form on a base

of mechanical analog of a viscoelastic material like a set

of rings and plungers known as Generalized Linear

Standard Solid (see [2]). Moreover, real data, both field

and laboratory, proved that (see [5-7]): A correct model-

ing scheme should yield a constant Quality Factor Q and

the corre sponding di spersion r e lation.

The common way to implement these properties of

real media in mathematical model is again just men-

tioned Generalized Standard Linear Solid (GSLS). For

this model the Hook’s law is written down as:

;

σσ

∑

iiR i

σε

σ τετ

∂∂



+=+



∂∂



and introduces Quality Factor as:

() ()

Lll

σε

εσ

ωτ τ

ωτ

ωωτ τ

ωτ

−+ +

=−

∑

The problem is how to choose a set of parameters

σε

ττ

providing desired behavior of Quality Factor

over a predefined interval of time frequencies.

We resolve it on the base of Least Squares techniques:

1 12

(,)|( )|minJQ Qd

σε ω

ττω ω

−−

= −→

∫

In order to find desired sets of relaxation times we ap-

plied

-method proposed in [1] and modified recently in

[3]. Its main advantage is essential reduction of a number

of relaxation times-for realistic values of Quality Factor

(

10Q>

) one can manage with single value of

(and

keep it the same through the target area!) and a couple o f

. Application of

-method for GSLS gives the fol-

lowing representation of Hook’s law:

(, )(1)(, )

1(,)

Rll

tx MLtx

Me xd

σ τε

τετ τ

−

=+−

−∑∫

Let us diffe rentiate this relation with respect to time

() ( )

(

)( )

lll

ML M

txex d

σσ

σττ

εετ τ

ττ

−

∂

∂=+−

∂∂



×−





∑∫

and introduce l-th memory variable

( )()

(, ),,]

rtxM txex

σσ

τε ετ

ττ

−

=−−



∫

(2)

Straightforward differentiation of (2) provides the fol-

lowing equation for this variable:

(, )

r rtx

σσ

ττ

∂∂

=−−

∂∂

and we come to the following system of first-order PDE

for viscoe l a s t ic wave pr o pa ga tion:

∂∂



(1 )

ML r

στ

∂∂

=++

∂∂

∑

∂∂



=−+



∂∂



As one can see, implementation of GSLS claims extra

independent “memor y” variables per relaxation mechan-

ism per stress component. In particular, for 3D models 

mechanisms claim

new variables besides three dis-

placements and six stresses. This leads to significant in-

crease of memory demands for simulation of waves’

propagation.

3. Parallel Implementation

Parallel Implementation is based on repr esentation of the

area of computations as superposition of disjoints sub-

domains touching one another along some contact sur-

face. Each of them is assigned to specific Processor Unit

(PU). Waves traveling in the model pass through differ-

ent subdomains, which require communication between

neighboring processors. Es sentially different geometry o f

the target area used in Sonic Log and Cross-hole Tomo-

graphy claims necessity of different approaches to Do-

main Decomposition (DD) as well. Let us consider both

of them.

3.1. Sonic Log

Target area for this statement is a circular cylinder

around the well with essential prevalence of its length in

comparison with width-typical length is about 10/12 me-

ters, while radius of this target area is no more th an 1/15

G. RESHETOVA ET AL.

meters. Taking this into account we slice the total 3D

model into a number of disc-like subdomains Ω. Finite

difference scheme assumes communication between

neighboring processors requiring them to exchange func-

tion values on the interfaces between elementary discs.

It should be noted that DD on the base of this rather

simple geometry easily provides possibility to guarantee

uniform load of Processor Units involved in computa-

tions. Another advantage of the chosen DD is in ex-

tremely small portions of data PU should interchange at

each time step and, so, extremely small waiting period

before computation on the next time step would be done.

The MPI (Message Passing Interface) library is applied

for arranging the above-mentioned send/receive proce-

dures and special efforts are paid in order to minimize

idle time of Processor Units due to the data exchange. In

order to provide this we start computations for each sub-

domain from its interior widening them towards inter-

faces and use non-blocking functions Isend and Ireceive

in order to arrange data exchange between neighboring

PU.

Special attention was paid to analysis of effectiveness

and scalability of this approach. This analysis was per-

formed by the series of numerical experiments performed

on the cluster HKC-160 (Siberian Supercomputer Center,

Novosibirsk) made of 80 computation modules (hp Inte-

grity rx1620 with two PU Intel Itanium 2, each of 1.6

Ghz, 3 Mb cache, 4 Gb RAM) connected via 24-port

commutator InfiniBand (10 Gbot, Cluster Interconnect).

Peak performance of the cluster is about 1 Tflop/s.

In order to estimate effectiveness of parallelization , fixed

computational area was decomposed on different quanti-

ty of subdomains, so simulation was performed for the

same target area, but with increasing number  of PU.

For each effectivenes s is found as

4* (4)

() * ()

time

eff nntime n

(3)

where

()timen

is computer time expended by  PU

for simulation. We start with

4n=

in order to provide

from the very beginning the same amount of data ex-

change between adjacent PU. On Figure 1(a)) one can

see effectiveness computed for two different lengths of

computational area: 12 meters (circles) and 24 meters

(rectangles). These results are completely predictable -

the less is load of PU the less is effectiveness. It is con-

firmed by behavior of both curves - for target area of 24

meters the load of single PU is twice as many as for tar-

get area of 12 meters and decrease of effectiveness is not

so sharp, but finally they become the same.

Now let us perform the series of numerical experi-

ments in opposite manner-we fix size of elementary

subdomains but increase their quantity proportionally to

quantity of PU. This result one can see on Figure 1(b))

Figure 1. Scalability by numerical experimentations: Up-

per-Effectiveness with fixed computational area; Down-

Total time with fixed load of PU.

and certain that the computational time does not change

under increase of quantity of PU.

So, we can conclude, that the key parameter is ratio of

amount of data being load to each of PU and amount of

data it should exchange with its neighbor. The higher is

this ratio the higher is effectiveness of parallelization. In

particular, if this ratio is fixed it does not matter how

many PU are involved in simulation.

3.2. Cross-Hole Tomography

Geometry of computational area used for cross-hole to-

mography is very different in comparison with previous

one-now it is parallelepiped with length (distance be-

tween wells) about 300/500 meters, approximately the

same depth (vertical size) and rather narrow width of ≈

100 meters. Next, contrary to Sonic Log, now local grid

refinement should be implemented for two different spa-

tial locations—around wells with sources and receivers,

so we need to pay special attention in order to provide

uniform load of PU. The easiest way to do this is to per-

form straightforward application of described above 1D

Domain Decomposition along cross-well direction, but it

will lead to necessity of huge data exchange between PU

G. RESHETOVA ET AL.

produced by wide contact surface between adjacent sub-

domains which leads to essential loss of parallelization

effectiveness. In order to minimize amount of data ex-

changed between PU we need to reduce surface area of

subdomains under fixed number of grid points per each

PU—that is to choose them as close to cubes as possible.

Under this approach each PU interchange data with 3/6

neighbors.

4. Numerical Experiments

To conclude the paper let us briefly consider simulation

results for a couple of realistic models performed on the

previously descri bed cluster HKC-160.

Sonic Log

The series of numerical experiments have been imple-

mented for a range of source frequencies, positions and

models of surrounding elastic media. For illustration let

us consider the model with well completion and vertical

crack presented on the Figures 2 and 3 and possesses the

following structure:

1) Vertical borehole with radius 0.1 m filled with a mud

with

=1500 m/s,

1000=



kg/m3, Quality Fac-

tor

65Q=

;

2) Steel tube encircling borehole; its width is equal to

0.01 m, wave propagation velocities

5600

m/s,

3270

m/s, =7830 kg/m3, Quality Factor

100Q=

;

3) The casing around steel tube; its width is equal to

0.04 m, its elastic parameters are the following:

= 4200 m/s,

= 2425 m/s,



= 2400 kg/m3,

Quality Factor

80Q=

;

4) Background-homogeneous viscoelastic layer with

wave propagation velocities

= 4989 m/s,

2605 m/s,



= 2400 kg/m3, Quality Factor

100Q=

;

5) Background-homogeneous viscoelastic layer with

wave propagation velocities

= 3208 m/s,

1604 m/s,



= 2400 kg/m3, Quality Factor

60Q=

;

Figure 2. Radial and azimuthal grid refinement around wel l

completion.

Figure 3. Upper—Vertical cross-section of the model; Mid-

dle—Horizontal cross-section of the model; Down—Zoomed

version of Middle figure.

6) Background-homogeneous viscoelastic layer with

wave propagation velocities

= 2650 m/s,

1219 m/s,



= 2400 kg/m3, Quality Factor

15Q=

;

7) Vertical crack started at r = 0.18 m and finished at r =

0.55 m, with 0.02 m in width ((see images b) and c))

filled with the same mud as within the borehole.

On the Figure 4 one can see snapshots stress rr

G. RESHETOVA ET AL.

Figure 4. Snapshots for

propagating through the plane

(, )0π

. Left column-the model without crack; Right col-

umn-difference between the calculations for the model with and without vertical crack.

during its propagation within borehole and adjacent rocks

with and without of vertical crack. There is a series of

snapshots for components through the plane

(0,π)

for axial source position with dominant frequency of 10

kHz located in

0.5z=

m. Two vertical red lines are

projections of interface between borehole and its com-

pletion with surrounding rocks, black horizontal lines

indicate interfaces of the layers and blue lines follow the

crack. The left column figures illustrate wave propaga-

tion for the model without of crack, while the right col-

umn figures represent the propagation of difference be-

tween the wave field with and without of crack.

5. Conclusion

The current version of developed software for simulation

of wave propagation in 3D heterogeneous viscoelastic

multiscale media opens possibility of careful study of a

rather wide range of realistic problems connected with

borehole-based geophysics. Of course, this version is not

perfect and should be further developed in different di-

rections. In particular, one of the most challenging prob-

lems is to perform full scale simulation for media with

extremely different scales, especially with microstructure

like caverns, pores and so on. In order to do that with

reasonable computer cost, one should be able to refine

G. RESHETOVA ET AL.

locally not spatial grid cells only, but time step used for

computations as well. Another important axis of devel-

opment is connected with necessity to provide an oppor-

tunity to perform simulation for more realistic mechani-

cal models and above all media with anisotropy induced

by fine layering, clusters of microcracks and stresses

around a well.

6. Acknowledgements

This work was supported by the Russian Foundation for

Basic Research (projects no. 11 -05-00947, 13-05-00076

and 13-05-12051).

REFERENCES

[1] J. O. Blanch, J. O. A. Robertsson and W. W. Symes ,

“Modeling of a Constant Q: Methodology and Algorithm

for an Efficient and Optimally Inexpensive Viscoelastic

Technique,” Geophysics, Vol. 60, No. 1, 2005, pp. 176-

184.

[2] R. M. Christensen, “Theory of Viscoelasticity—An In-

troduction,” Academic Press, Waltham, 1982.

[3] S. Hestholm, et al., “Quick and Accurate Q Parameteriza-

tion in Viscoelastic Wave Modeling,” Geophysics, Vol.

71, No. 5, 2006, pp. T147-T150.

[4] V. I. Kost in, D. V. Pissarenko, G. V. Reshetova and V. A.

Tcheverda, “Numerical Simulation of 3D Acoustic Log-

ging,” Lecture Notes in Computer Sciences, Vol. 4699,

2007, pp. 1045-1054.

http://dx.doi.org/10.1007/978-3-540-75755-9_122

[5] F. J. McDonal, F. A. Angona, R. L. Mills, R. L. Sengbush,

R. G. van Nostra nd and J. E. White, “Attenuation of Shear

and Compressional Waves in Pierre Shale,” Geophysics,

Vol. 23, No. 3, 1958, pp. 421-439.

http://dx.doi.org/10.1190/1.1438489

[6] S. A. Siddiqui, “Dispersion Analysis of Seismic Data,”

M.S. Thesis, University of Tulsa, Tulsa, 1971.

[7] P.C. Wuenschel, “Dispersive Body Waves—An Experi-

mental Study,” Geophysics, Vol. 30, No. 4, 1965, pp.

539-551. http://dx.doi.org/10.1190/1.1439620