PROCEEDINGS of the23rd International Congress on Acoustics
9 to 13 September 2019 in Aachen, Germany

Parallel implementation for three-dimensional acoustic field computation
in a penetrable wedge by image source methodWenbin XIAO 1; Yongxian WANG 1*; Wei LIU 1; Qiang LAN 1; Xinghua CHENG 1; Zijie ZHU 1;
Xin WANG 1; Ben LUO 2, Dezhi WANG 1; Jiani WU 1, Lilun ZHANG 1
1

College of Meteorology and Oceanology, National University of Defense and Technology, P.R.China
2

School of Traffic & Transportation Engineering, Central South University, P.R.China

ABSTRACT
Present work aims to the parallel implementation of message passing interface (MPI) and open
multi-processing (OpenMP) for the three-dimensional acoustic field computation in a penetrable wedge.
Possible practical application of the wedge waveguide propagation is extended to the continental shoulder orshallow water environment. By the image source method, the acoustical propagation solution is inconjunction with an extension of reflection analysis of spherical wave filed from a planar interface. Multiple
reflections at the boundaries are taken as the originating from an image source depended on the surface andbottom positions. Analytical expression of acoustical pressure here is coupled with the Bessel function
expansion in complex form. Two range-dependent cases, such as the along the slope (two dimension) andcross the slope (three dimension), are tested by the MPI and OpenMP parallelization and performance
optimization. The multi-core platform of “Knights Landing” or KNL and many-core platform of “Tianhe-2”
are applied to the utilization. The performance gain observed in using OpenMP of the FORTRAN code onhigh performance computing facilities is demonstrated.
Keywords: wedge, underwater acoustic propagation, MPI, OpenMP, Tianhe-2, KNL

1. INTRODUCTION
Underwater acoustic propagation is one of the basic contents in computational ocean acoustics [1]. Earlystudies mainly focused on the two-dimensional problem for the sake of limitation by the ocean environment
data and computing ability. However, the underwater acoustic effects such as interface reflection greatly
influence the propagation in a real complicated environment of ocean. In the process of solving thethree-dimensional underwater acoustic field, some key difficulties are demonstrated as the highcomputational complexity and long time consuming, which limits its application in the matched field
location, geoacoustic property inversion and array performance assessment for the calculating speed
requirement [2]. With the rapid development of computer hardware and parallel optimization technology, it
is one of the effective methods to propose high-efficiency and large-scale parallel algorithm of
three-dimensional underwater acoustic propagation based on the latest high-performance platforms.
The three-dimensional acoustic propagation in a penetrable wedge is a theoretical model extracted fromthe shallow continental shelf and slope oceans [3]. An analytical solution of the sound field in a penetrable
wedge with a stratified fluid or elastic basement was proposed by Dean and Buckingham [4], in which thesurface and bottom reflections are explored by the image source method. Luo et al presented a generalized
coupled-mode formulation was presented to calculate the wedge acoustic field by applying direct global
matrix approach to obtain the modal expansion coefficients [5]. By using a square-root Helmholtz operator
splitting algorithm, Ying et al proposed an alternating direction implicit three-dimensional fluid parabolic
equation solution method with enhanced accuracy for the wedge problem [6]. The adiabatic mode parabolic
equations were also used to approximate the solution of the three-dimensional Helmholtz equation bymodal decomposition of the acoustic pressure field, in which the mode amplitudes satisfy parabolic
equations that admit analytical solutions in the special case of the three-dimensional wedge [7]. Tang et alobtained the characteristics of horizontal refraction in the wedge-shaped ocean and released the MATLAB*1Correspondence: yxwang@nudt.edu.cn

1921code on the website of Ocean Acoustics Library [8]. Present work is aimed at the parallel implementation of
three-dimensional acoustic field in a wedge waveguide after the transformation Tang’s code to theFORTRAN version and performance test on the “Tianhe-2” and KNL platforms. It is expected to support
the benchmark of three-dimensional range-dependent propagation problem from the perspective ofcomputational speed.

2. BASIC MATHEMATICL PROBLEM
2.1

Reflection of Plane Waves
If the interface between two media is untitled as Figure 1, the spherical wave field is decomposed into alinear superposition of plane waves. For the sake of boundary reflection, these plane waves are only
modified with amplitude, phase and direction of propagation. The reflection coefficient V on the interface
is introduced to account for the quantitative change of plane wave properties. In the local CartesianሬԦ is a unit vector normal to the interface. The wave-number vector Ɉcoordinate system of image source, ‫܍‬ሬԦሬԦ are here intersected at the incidence angle ߮, which can be described as below.
and ‫܍‬ܿ‫ ߮ ݏ݋‬ൌԦ ‫܍ڄ‬ሬԦ‫ܓ‬(1)

௞In frequency domain, the spherical wave field at the receiver point is the contribution of homogeneousand inhomogeneous plane waves. For the former plane wave, the incidence angle ߮ is real from Ͳ toߨȀʹ. For the latter one, ߮ is varying form ߨȀʹ to ߨȀʹ െ ݅∞ shown in Figure 2. Taken into account ofplane-wave reflection, the coefficient V can be demonstrated as following.
ܸሺ߮ሻ ൌሺఘమ Ȁఘభ ሻఊభ ିఊమ(2)ሺఘమ Ȁఘభ ሻఊభ ାఊమ

Figure 1 – Image source for untitled interface

Figure 2 – Range value of ߮ for untitled interface

For the formula above, ߛଵ and ߛଶ are the magnitudes of the wave-number vector normal to the media
interface in the upper and lower fluids respectively, whose expressions are responding to thatߛଵ ൌఠܿ‫߮ ݏ݋‬,௖భఠ௖మభଶߛଶ ൌ ௖ ටͳ െ ቀ௖మ ቁ ሺͳ െ ܿ‫ ݏ݋‬ଶ ߮ሻ(3)where ߩଵ and ߩଶare the density values in the upper and lower fluids respectively, ܿଵ and ܿଶ are the soundspeeds in the upper and lower fluids respectively and ߱ is angular frequency for the harmonic point
source.
2.2

Image Source for an Inclined Interface
In the wedge-like waveguide problem, there is an angle ȣfor the interface between two media, which
can be drawn as the inclined case as Figure 3. In order to decompose the plane waves, it is necessary toapply the spherical polar coordinate system for the spherical wave field. Then the wave-number vector ɈሬԦand coordinates ሺ‫ݔ‬௥ ǡ ‫ݕ‬௥ ǡ ‫ݖ‬௥ ሻ at receiver point can be transformed into the spherical polar forms.
݇௫ ൌ ݇ ‫߶ ݏ݋ܿ ߠ ݊݅ݏ‬, ݇௬ ൌ ݇ ‫߶ ݊݅ݏ ߠ ݊݅ݏ‬, ݇௭ ൌ ݇ ܿ‫ߠ ݏ݋‬(4a)‫ݔ‬௥ ൌ ܴ ‫ߦ ݏ݋ܿ ߞ ݊݅ݏ‬, ‫ݕ‬௥ ൌ ܴ ‫ߦ ݊݅ݏ ߞ ݊݅ݏ‬, ‫ݖ‬௥ ൌ ܴ ܿ‫ߞ ݏ݋‬(4b)
ሬሬԦ and the ‫ ݖ݋‬axis and ߶ or ߦ is the angle between theሬԦ or ‫܀‬
where ߠ or ߞ is the angle between ɈሬሬԦ onto the ‫ ݕ݋ݔ‬plane and the ‫ ݔ݋‬axis.
projection of ሬɈԦ or ‫܀‬

1922

Figure 3 – Image source for inclined interface
By the coordinate transformation and variables separating, the mathematical representation of reflected
acoustical field at the receiver point can then be written in the following formula.
௜௞గȀଶି௜ಮ

‫݌‬௥ ൌ ଶగ ‫׬‬଴ଶగ݁ ௜ஐ೥ ሺఏሻ ‫ ߠ ݊݅ݏ‬ൈ ‫׬‬଴ ܸሺ߮ሻ݁ ௜ஐሺఏሻ ௖௢௦ሺథିకሻ ݀߶݀ߠ(5)ሬሬԦ.
where ȳ௭ ሺߠሻ ൌ ܴ݇ ܿ‫ ߞ ݏ݋ܿ ߠ ݏ݋‬and ȳሺߠሻ ൌ ܴ݇ ‫ ߞ ݊݅ݏ ߠ ݊݅ݏ‬are derived from the dot product of ɈሬԦ and ‫܀‬And the incidence angle ߮ is then expressed as the result of spherical components.
ܿ‫ ߮ ݏ݋‬ൌ ‫ ݊݅ݏ ߠ ݊݅ݏ‬ȣ ܿ‫ ߶ ݏ݋‬൅ ܿ‫ ݏ݋ܿ ߠ ݏ݋‬ȣ
(6)
Here ߮ is the function of both ߠ and ߶, which adds the integral difficulty of formula (5). It is wellknown that the integration of exponential term can be expanded in a Bessel function series. In theassumption of isotropic media, the reflected acoustical pressure can be simplified as below in the chosencoordinate system.
‫݌‬௥ ൌ௜௞ గȀଶି௜ಮ ௜ஐ ሺఏሻ௔ ௃ ሺஐሻஶఔఔ݁ ೥ ‫ ߠ ݊݅ݏ‬ቀ బ ଶబ ൅ σஶ‫׬‬ఔୀଵሺെͳሻ ܽଶఔ ‫ܬ‬ଶఔ ሺȳሻ ൅ ݅ σఔୀ଴ሺെͳሻ ܽଶఔାଵ ‫ܬ‬ଶఔାଵ ሺȳሻቁ ݀ߠଶగ ଴(7)గܽఔ ൌ ʹ ܿ‫ݏ݋‬ሺ ߥߦሻ ‫׬‬଴ ܸሺ߮ሻ ܿ‫ݏ݋‬ሺ ߥ߶ሻ݀߶(8)where ‫ܬ‬଴, ‫ܬ‬ଶఔ and ‫ܬ‬ଶఔାଵ are the first-kind Bessel functions of zero order, ʹߥ order and ሺʹߥ ൅ ͳሻ orderrespectively and ܽఔ is the Fourier coefficient.
2.3

Acoustic Field in a Three-dimensional Wedge
The geometry of a wedge-like waveguide in three-dimension acoustical problem is shown as Figure 4.
The angle of wedge apex is defined as ȣwith the sea surface and bottom. A harmonic point source isplaced in the three-dimension acoustical filed. In an orthogonal coordinate system, the origin ܱ is fixed at
the wedge apex, the positive ܱܺ axis is toward the wedge volume, the positive ܱܼ axis is pointing
downward and ܱܻ axis is along the wedge slope. Present work is aimed at the numerical calculation of
acoustical propagation at receiver points in the three-dimension wedge-like waveguide.

Figure 4 – Geometry of a wedge-like waveguide
Each reflection on both the surface and bottom boundaries raises the number of image source, which islocated on the circle centered at the wedge apex. There are totally four kinds of image sources at specificreflection number ݊௕ from the bottom boundary. In order to identify the image source uniquely, aparameter ݈ is adopted with values of ሺെʹǡ െͳǡͳǡʹሻ. Whether the first reflection occurred at the surface orbottom is related to ݊௦ , which is the number from the surface reflections. In the case of ݈ ൏ Ͳ, the image

1923

source is below the bottom. For ݈ ൐ Ͳ, the image source is above the bottom. The angle between a line
passing through the source and wedge apex is defined as ߙ௦ ൌ ȣ െ ‫ି݊ܽݐ‬ଵሺ ‫ݖ‬௦ Ȁ‫ݔ‬௦ ሻ. If both the source
location and wedge angle are specific, the angle of image source can be denoted as Φ௡௟ , which is taken theinterface-reflection effect into account.

) nl­ D s  24( nb  1),°°D s  24nb ,®
°D s  24nb ,° D  24( n  1),b
¯ sl2l 1l

1l

2(9)

The geometry of the image source in acoustic field is shown in Figure 5. For discrete receiver points inthe space, the complex pressure of acoustic field can be obtained by the double integration as formula (10).
If only the actual source is taken into account, the acoustic pressure at receiver point is related to the direct
wave. However, the reflection of plane wave on the surface or bottom interface exerts significant effect on
the acoustic field. For certain bottom reflection, more exact pressure should be the superposition by fourkinds of image source. All the image sources are distributed on the arc such as Figure5. The quantitativevalue of acoustic field by a large enough number݊௕ is less marked with the increase of bottom reflection,
which hints that the acoustic pressure is convergent in the computational field.
௡σሬሬԦǡ Φ௡௟ ሻσ௟ ‫݌‬௥ ሺ‫܀‬‫݌‬ሺ‫ݔ‬௥ ǡ ‫ݕ‬௥ ǡ ‫ݖ‬௥ ሻ ൌ σ௡೘ೌೣ(10)್ ୀ଴Figure 5 – Images of the source in acoustic field

3. DESIGN OF PARALLEL STRATEGY
3.1

Analysis of Computational Flow
Computational flow of the acoustic field by image source method is illustrated in Figure 6. Firstly,
three-type parameters involving the actual source, receivers and environment are input. For actual source, itis necessary to obtain the frequency and coordinates in space. The exact position of receiver points shouldbe also specific. Parameters of the environment are inclusive of the angle of wedge apex, sound speed inwater and bottom.
Three loops are included in the whole framework. The number of bottom reflection݊௕ is related to theouter loop, which initiates the calculation or update of coefficient V for plane-wave reflection as part one.
For the reason that specific image source is corresponding to the angle Ȱ௡௟ , the middle loop is related tothe different kind of image source. If the interface is inclined, the Cartesian coordinates of actual source orimage source and receiver points are then transformed into the cylinder ones. The output of coordinatetransformation is then input to calculate the first-kind Bessel functions with multiple sequences denoted aspart two. The contribution of actual source or image source to the acoustic field is described in aquantitative way by the production of Fourier coefficient and Bessel function in part three, where thereflection coefficients are needed from the part one of outer loop. Therefore, the computation of inner loop
depends on the maximum number of Bessel sequence. Based on the superposition principle in numericalintegration, the acoustic pressure is obtained by the Simpson integration for integral variableߠ, which ismarked as part four. If the difference between the former pressure at ሺ݊௠௔௫ െ ͳሻ and present one at ݊௠௔௫is small enough, the numerical results can be considered to be convergent and the whole computation is1924

then ended. From the view of computational complexity, the four parts mentioned above are very
time-consuming.

Figure 6 – Computational flow of the acoustic field by image source method
3.2

Task-level Decomposition
From the analysis of computational flow above, the acoustic pressure is solved by iteration related to
the number of reflection, which determines that the outer loop can’t be parallel implemented. Additionally,
the reflection coefficient is relevant between the part one and part three. A feasible parallel solution may bemined for the fusion of middle and inner loops, in which the degree of parallelism is only the production
between the maximum number of Bessel sequence and four kinds of image source. Such parallel
implementation is limited to the high scalability.
The task-level decomposition in coarse granularity is here explored from the receiver points illustrated
in Figure 7. In the condition of specific parameters for the actual source and environment, the calculation ofacoustic pressure at each receiver point is absolutely independent. Therefore, the space of receivers can bedivided into many subspaces in three dimensions and each subspace is then mapping to the specificcomputational task. This parallel implementation can be easily decomposed in the task level and extended
to many CPU cores by the MPI model.

Figure 7 – Mapping from the subspaces to the computational tasks
The calculation of large matrix in some parts is especially compute-intensive, which may besuitable for the parallel implementation of OpenMP model. With the part one as example, pseudo codeby the OpenMP parallel is illustrated in Figure 8. The dimension of reflection coefficient matrix is ߠ͖߶, which are corresponding to the integral variables in formula (7). For the sake of computationalaccuracy, the discrete number of ߠ or ߶ is generally set in the thousands. From the formulas (3) and(6), the variable ߮ is the function of ߠ and ߶ and each element in the reflection coefficient matrix is
relatively independent, which is applicable to the task decomposition in different threads. During theprogram execution, each subtask is mapping to the computation of specific thread according to the setof OpenMP environment variable.

1925

!$omp parallel &!$omp & default(shared) &!$omp & private( iph, ith, cos_phi_b, gam1, gam2)!$omp do collapse(2)do iph = 1, nph% nph is the discrete number of ߶
do ith = 1, nth% nth is the discrete number of ߠ
cos_phi_b = sin(theta_b)×sin(theta(ith))×cos( phi(iph))+cos(theta_b)×cos( theta(ith)) % formula (6)gam1 = k×cos_phi_b
gam2 = k×sqrt( (c1/c2)2-1.0 + cos_phi_b2 )% formula (3)BotRefl = (rho2_rho×gam1-gam2)/(rho2_rho1×gam1+gam2)% formula (2)Vtotal(ith, iph) = Vtotal(ith, iph)×BotRefl% update V
enddo
enddo
!$omp end do!$omp end parallel

Figure 8 – Pseudo code of the OpenMP implementation for reflection coefficient

4. PERFORMANCE OF PARALLEL OPTIMIZATION
4.1

Platforms and Cases
Two types of computer platforms are considered in present work. They are the many-core platform of“Tianhe-2” and multi-core platform of Intel Xeon Phi 7250 (formerly Knights Landing, KNL). Super
computer “Tianhe-2” is typical cluster structure, whose computing nodes are connected by high-speed
internet and can communicate with each other. It is configured with the tens of thousands of Ivy Bridge
processor and Xeon Phi coprocessor. This type of test platform was ranked as the fastest supercomputer in
Top 500 from 2013 to 2016 [8]. KNL is the second generation with Many Integrated Core architecture,
which is specifically designed for highly parallel workloads by Intel [9]. The working mode of KNL may beused as the master processor or a coprocessor whose peak performance per node can reach 3.0T FLOPSwith double precision. KNL-7250 is advantageous of the optimal performance to power ratio between thefour products launched by Intel in 2016. Table 1 lists the hardware and software configurations of above
two platforms at one node.
Table 1 – Configurations of the test platforms at one node
Items

Tianhe-2

KNL

CPU typeCPU cores
CPU frequency
Main memory
Operating system

Intel Xeon E5-2692 v2
242.2 GHz64 GB
Red Hat 6.5

Intel Xeon Phi 7250
68
1.4GHz
16 GB
Red Hat 4.8

The test cases of three-dimensional acoustic filed in a penetrable wedge include two conditions of thealong-slope and cross-slope propagations at certain water depth. That is to say, the first case is for thereceiver points on the ܱܼܺ plane and the second one is parallel to the ܱܻܼ plane in Figure 4. Some
identical settings in both cases are inclusive of source frequency݂ ൌ ʹͷ , wedge angle ȣ ൌ ʹǤͺ͸ʹͶι,sound speed in water ܿଵ ൌ ͳͷͲͲȀ, sound speed in bottom ܿଶ ൌ ͳ͹ͲͲȀ, water density ߩଵ ൌͳǤͲȀ ଷ, bottom density ߩଶ ൌ ͳǤͷȀ ଷ and bottom attenuation ͲǤͷȀߣ. Other parameter settings are
all listed in Table 2.
Table 2 – Settings of some parameters in the two test cases

1926

Parameters

First Case

Second Case

Receiver position
Discrete number of receiver points
Source position
Discrete number for integral ߠDiscrete number for integral ߮
Maximum number of bottom reflection
Sequence number of first-kind Bessel function
Number of image sourceሺ‫ݔ‬௥ ǡ Ͳǡ ͵Ͳሻ512ሺͶͲͲͲǡ Ͳǡ ͳͲͲሻሺͶͲͲͲǡ ‫ݕ‬௥ ǡ ͵Ͳሻ512ሺͶͲͲͲǡ Ͳǡ ͳͲͲሻ70021001
92387002
1001
3050122

4.2Parallelization of Sequential Code
It is very important to parallelize sequential code for that the code must be operated in sequence in theparallel zone, single process or single thread. In order to compare the performance by different compiler
and its options, several schedules are explored for the first case of wedge acoustic field on Tianhe-2 andKNL, whose results are illustrated in Figure 9. The high-level optimization of compiler is inclusive of –O3and –O2 and the low-level optimization is for –O1 and –O0. The larger value of –O means higher
parallelization. The option of –funroll-all-loops is loop unrolling, in which the compiler unrolls the loops toreduce the number of branches. Test results show that the consuming time on Tianhe-2 is rather less thanthat on KNL. The parallel performance is better if the compiler options of –O3 and –funroll-all-loops areboth combined for the ifort than gcc.

Figure 9 – Parallelization of sequential code by different compiler and options
4.3Multi-thread Parallelization
Each node of KNL has 68 CPU(s) which can launch four threads simultaneously. For the first andsecond cases, multi-thread parallelization is implemented on the KNL platform, whose results are shown inFigure 10. If the thread number increases to 64, the speed up of performance is excellent. However, theparallelization is not good if the thread number is larger than 64. At these conditions, the multi-threadparallelization is up to the load limit of KNL, which is corresponding to the half-full-load status. However,the maximum number of threads is 272 for KNL in physical.

(a) First case(b) Second caseFigure 10 – Multi-thread parallelization on KNL platform
4.4Hybrid MPI/OpenMP Parallelization
The performance test of hybrid MPI/OpenMP is also explored for the two cases. On the Tianhe-2
platform, different combination is configured for the 1 h24, 2 h12, 4h6, and 8h 3, where the former
number is corresponding to the processor and the latter one is to thread. From the test results inFigure 11, the performance of large thread number is better than that of processor number.
Additionally, the calculating task of part three is very significant for the second case. By the hybrid
MPI/OpenMP test, it can improve the performance by the proper combination of process and threadnumbers.

1927

(a) First case(b) Second case
Figure 11 – Hybrid MPI/OpenMP parallelization on Tianhe-2 platform

5. CONCLUSIONSFor the three-dimensional acoustic field computation in a penetrable wedge by image source
method, the parallel model by MPI and OpenMP is proposed in present work. The time-consuming andspeed-up performances of sequential parallel, multi thread and hybrid MPI/OpenMP parallels are also
tested on the Tianhe-2 and KNL. Some results demonstrate the good acceleration and scalability forfurther exploration.

ACKNOWLEDGEMENTS
Present work received financial support from the National Natural Science Foundation of China(No.51709267, No.61702531, No.61379056, No.61502516), National Key R & D Program of China(No.2016YFC1401800), Natural Science Foundation of Hunan Province, China (No. 2017JJ2305) andthe Scientific Research Project of NUDT (No.ZK17-03-31, No.ZK16-03-31, No.ZK16-03-46).

REFERENCES
1. Jensen F B, Kuperman W A, Porter M B, et al. Computational Ocean Acoustics (Second Edition).
Springer New York Dordrecht Heidelberg London, 2011.
2. Xu Min, Wang Yongxian, Anthony T C, Yue Hao. Parallelization and performance optimization ofcalculation in three-dimensional underwater acoustic propagation on modern many-core processor. IEEEInternational Conference on Signal Processing, Communications and Computing, 2018.
3. Kevin D H, Richard L C, James J M. Comparison of hybrid three-dimensional modeling withmeasurements on the continental shelf. Journal of the Acoustical Society of America, 2012. p.
1680-1688.
4. Deane G B, Buckingham M J. An analysis of the three-dimensional sound field in a penetrable wedge
with a stratified fluid or elastic basement. Journal of Acoustical Society of America, 1993. p. 1319-1328.
5. Luo Wenyu, Yang Chunmei, Zhang Renhe. Generalized coupled-mode formulation for soundpropagation in range-dependent waveguides. Chinese Physics Letters, 2012. p. 14302.
6. Ying Tsong Lin, Jon M Collis, Timothy F Duda. A three-dimensional parabolic equation model of soundpropagation using higher-order operator splitting and Pade approximants. Journal of the AcousticalSociety of America, 2012.p. 364-370.
7. Pavel S Perov, Frederic Sturm. An explicit analytical solution for sound propagation in athree-dimensional penetrable wedge with small apex angle. Journal of the Acoustical Society of America,
2016.p. 1343-1352.
8. Tang Jun, Piao Shengchun, Yang Xianpeng, Zhang Minghui. An analysis of horizontal refraction in theproblem of sound propagation in wedge-shaped oceans. International Conference of Oceans, 2017.
9. Tianhe-2 championed global Top500 for sixth time. Sino Cast Daily Computers & Electronics Beat,2015.
10. Jim Jeffers, James Reinders. Intel® Xeon Phi™ Coprocessors High-Performance Programming Knights
Landing Edition. San Francisco: Margan Kaufmann, 2016.

1928