PROCEEDINGS of the
23rd International Congress on Acoustics
9 to 13 September 2019 in Aachen, Germany

6JTGGFKOGPUKQPCNCEQWUVKERCTCDQNKEGSWCVKQPOQFGNDCUGFQP)27
RTQEGUUKPIKeunhwa Lee1; Woojae Seong2; Youngnam Na312
3Sejong University, South Korea

Seoul National University, South Korea

Agency for the Defense Development, South Korea

ABSTRACT
We introduce a three-dimensional acoustic parabolic equation (PE) model based on GPU processing. Thismodel is developed by splitting the full exponential operator in the three-dimensional split-step marchingsolution in the Cartesian coordinate. Each exponential operator is approximated by the Taylor series and thesum representation of Padé approximation, which gives an opportunity for GPU parallel processing of the 3DPE solver. An advantage of this model is to implement higher-order cross-terms appearing to the operator
approximation efficiently. In the talk, we will show the computational performance of GPU based PE codewith the comparison of a traditional code for several ocean environments
Keywords: Three-dimensional parabolic acoustic equation (3DPE), GPU, Parallel algorithm

1. INTRODUCTIONRecently, some of the three-dimensional acoustic parabolic equation (3D PE) models considering
the cross-terms in the two-dimensional (2D) square root operator have been presented. While thesemodels have advantages in the accuracy and the ܿ଴ -insensitivity compared with the existing 3D PE
model without the cross-terms, the total calculation time increases more due to the calculation of thecross-terms and the use of filter. Fortunately, the PE algorithm is known to be compatible with theparallel processing inherently.
In this proceeding, in order to improve the efficiency of the 3DPE model, we present a 3DPE
formulation proper to the GPU processing and the computational strategy. Specifically, we introduce
the 3DPE model based on Lee et al. approach and its parallelized formulation, known as the sum
representation of the operator approximation. Moreover, we suggest a strategy for the matrix
multiplication and inversion using the cuSparse library in CUDA/C. With these two approaches, thetotal calculation performance increases about 5 ~ 10 times.

2. METHOD
2.1 3D PE solver with cross-terms in the Cartesian coordinateWhen the marching direction is set to be x-axis, the 3D PE solver with cross-terms is given byp
1  Di , pY p 1  Di , p Z pyz(1)u ( x  'x) ¦ i u( x) .
i 1 1  E i , pY i 1 1  E i , p Z i 0Here, the variable u is defined by u
1nasalkh2@sejong.ac.kr
wseong@snu.ac.kr
3ynna@add.re.kr
21910

pe jk0 x / D with the acoustic pressure p, the referencewavenumber ݇଴, and the energy-conserving parameter Ƚ.  is the order of Padé approximation,ߙ௜ǡ௣ and ߚ௜ǡ௣ are the Padé coefficients and ‫݌‬௬௭ is the order of Taylor approximation. i denotes theith order cross-term of Taylor series. Here, Y and Z operators are defined asU w §1 w ·(2)Y k02D¸,¨
D wy © U wy ¹

U w §1 w ·(3)D ,D wz ¨© U wz ¸¹
where n k ( y, z) / k0 is the refractive index where k is the complex wavenumber. These operators
can be transformed into the tridiagonal matrix using the second-order finite difference scheme.
Figure 1 illustrates the marching algorithm of the 3D PE solver, regarded as the solution of theinitial value problem. When the pressure field at  ൌ  is known, the pressure field at the next step
can be predicted with Eq. (1). In Eq. (1), the multiplication operator results from the Padé
approximation of the exponential operator to the product representation form.
The sum representation of Eq. (1) is derived asnn§
·Dˆ Y ·§Dˆ Z · § nyz nyzu ( x  'x) ¨1  ¦ i , n ¸¨1  ¦ i , n¸ ¨¨1  ¦¦ cˆij , nY i Z j ¸¸ u ( x) ,(4)ˆˆ¨¸¨¸i 1 1  E i , nY ¹©
i 1 1  Ei , n Z ¹ ©
i 1 j 1¹©
where  is the order of Padé approximation in the sum representation, ߙො௜ǡ௣ and ߚመ௜ǡ௣ are the Padé
coefficients and ݊௬௭ is the order of Taylor approximation. ܿƸ௜௝ǡ௡ is the coefficient of Taylor series.
In figure 2, the flow chart for the calculation of Eq. (4) is given. For convenience, the precise flowchart for the cross-term operation is omitted. Operations in each bracket of Eq. (4) can be performedin parallel. Theoretically, the calculation time in each bracket may be reduced by approximatelyଶ. This kind of parallelization concept has been already applied into several 2D PE orͳȀor ͳȀ݊௬௭
3D PE solvers.
Z(n2  1)  k02

Figure 1 - Illustration of the 3D PE algorithm in the Cartesian coordinate.

Figure 2 – Flow chart of the 3D PE algorithm in the sum representation.

1911

Figure 3 - Grouping of rows and columns for parallel processing.

More parallelization can be achieved for each tri-diagonal matrix multiplication and inversion infigure 2. For example, for the calculation of (1  Eˆi , nY )1 (Dˆi , nY )u( x) , the nz number of tri-diagonal
matrix multiplication and inversion are required when the size of the discretized ሺሻ is set to benz u ny . Since the matrix calculation for each row of ሺሻ is independent one another, theseoperations can be parallelized like those for Padé order. The same is true for the Z operator. As
shown in figure 3, we assume to divide the u(x) matrix into the sub-matrices. When the numbers of
the sub-matrix are set to be ܰ௦௭ in the z-direction and ܰ௦௬ in the y-direction, totally the ሺ݊ ൈ ܰ௦௭ ሻor ሺ݊ ൈ ܰ௦௬ ሻ loops will be performed in parallel. Moreover, the inversion of a tri-diagonal matrix
can be fast implemented with parallel algorithms like parallel cyclic reduction and recursivedoubling.
The CUDA/C library is used to code the suggested 3DPE algorithm in parallel. The matrix
inversion is performed by the function of cuSparse<t>gtsvStrideBatch in CUDA Toolkit 4.2
cuSparse library. The solver runs on a desktop with Intel i7-6700k quad-core processor, 56 GB RAM
Intel and a GeForce GTX 1080 GPU with 20 multiprocessors.

3. CONCLUSIONSIn this proceeding, the parallelized 3DPE solver using GPU based on CUDA library is presented.
The parallelization is carried out through two steps-the sum representation of the operator
approximation and the grouping of vectors in the matrix. With the suggested algorithm, the execution
time of the 3DPE is shortened to 1/5 ~ 1/10.

ACKNOWLEDGEMENTS
This work was supported by the Agency for Defense Development, South Korea, under Grant No.
UD170018DD, and by the KRISO, South Korea.

REFERENCES
1. Lin Y, Collis J. M, Duda T. F. A three-dimensional parabolic equation model of sound propagation using
higher-order
operator splitting
and Pade approximationts.
JAcoustSoc Am.
2012;132(5):EL364-EL370.
2. Sturm F. Leading-order cross term correction of three-dimensional parabolic equation models. J Acoust
Soc Am. 2016;139(1):263-270.
3. Castor K, Sturm F. Investigation of 3D acoustic effects using a multiprocessing parabolic equation basedalgorithm. J. Comput. Acoust. 2008;16(2):137-162.
4. Lee K, Seong W, Na Youngnam. Three-dimensional Cartesian parabolic equation model with higher
order cross-terms using operator splitting, rational filtering, and split-step Padé algorithm. J Acoust Soc
Am. 2019 (accepted).
5. NVIDIA, CUDA Toolkit 4.2 CUSPARSE Library. PG-05329-041_v01. February 2012.

1912