- Email: [email protected]
Parallel Acoustic Field Simulation with Respect to Scattering of Sound on Local Inhomogeneities
Available online at www.sciencedirect.com
ScienceDirect
This space is reserved for the Procedia header, do not use it This space is reserved for the Procedia header, do not use it Procedia Computer Science 108C (2017) 1753–1762
This space is reserved for the Procedia header, do not use it
International Conference on Computational Science, ICCS 2017, 12-14 June 2017, Zurich, Switzerland
Parallel Parallel Acoustic Acoustic Field Field Simulation Simulation with with Respect Respect to to Scattering of Sound on Local Inhomogeneities Parallel Acoustic Field on Simulation with Respect to Scattering of Sound Local Inhomogeneities 1 1 of Sound on Local Inhomogeneities AndreyScattering A. Chusov11 , Lubov G. Statsenko 1 , Alexsey P. Lysenko1 , Sergey N.
Andrey1 ,A. Chusov , Lubov G.1 ,Statsenko , Alexsey P. Lysenko , Sergey N.2 1 Kuligin 1 Nina A. Cherkassova 1 Petr P. Unru 1 , and Maya V. Bernavskaya 2 1 1 1 Nina A. Cherkassova Petr P. Unru , and P. Maya V. Bernavskaya Kuligin Andrey ,A. Chusov , Lubov G. ,Statsenko , Alexsey Lysenko , Sergey N. 1 1 1 2 Far Eastern Federal University, Vladivostok, Russia Kuligin1 , Nina A. 1 Cherkassova , Petr P. Unru , and Maya V. Bernavskaya Far Eastern Federal University, Vladivostok, Russia 2 2 2
[email protected] [email protected] 1 Saint-Petersburg Polytechnic University, Saint-Petersburg, Far Eastern Federal University, Vladivostok, Russia Russia Saint-Petersburg Polytechnic University, Saint-Petersburg, Russia [email protected] [email protected] [email protected] Saint-Petersburg Polytechnic University, Saint-Petersburg, Russia [email protected]
Abstract Abstract The report presents a developed approach to simulation of acoustic fields in enclosed media. The report a developed to simulation of acoustic fields enclosed sources media. Abstract This methodpresents is based on the useapproach of Rayleigh’s integral for calculation of in secondary This method is based on the use of Rayleigh’s integral for calculation of secondary sources generated a wavea falling ontoapproach media boundaries. The is media. highly The reportby presents developed to simulation of implementing acoustic fields algorithm in enclosed generated by aimplies wave falling onto media boundaries. The implementing algorithm is highly parallelizable, loosely coupled parallel branches with only few points of inter-thread This method is based on the use of Rayleigh’s integral for calculation of secondary sources parallelizable, implies loosely coupled parallel branches with only few points of inter-thread communication. On the other hand, the algorithm isThe exponential upon an average generated by a wave falling onto media boundaries. implementing algorithm isnumber highly communication. On the other hand, the algorithm is exponential upon an average number of reflections which occur to a single wave element emitted by a primary source, although parallelizable, implies loosely coupled parallel branches with only few points of inter-thread of reflections which occurthis to a single can wave element emitted byprovide a primary source, although for practical applications reducedis enough to results with communication. On the othernumber hand, the be algorithm exponential uponaccurate an average number for practical applications this number can be reduced enough to provide accurate results with reasonable time and space of reflections which occur consumptions. to a single wave element emitted by a primary source, although reasonable time and space consumptions. for The practical applications thisis number reduced enough to provideofaccurate results proposed algorithm based oncan thebeapproximate superposition acoustical fieldswith and The proposed algorithm is based on the approximate superposition of acoustical and reasonable time and space consumptions. provides adequate results, as long as the used equations of acoustics are linear. To fields calculate provides adequate results, as long as the used equations of acoustics are linear. To calculate scattering properties of reflecting boundaries, the algorithm representsofa acoustical geometric fields modeland of The proposed algorithm is based on the approximate superposition scattering properties of reflecting boundaries, the algorithm represents a geometric model of sound media propagation set of flat vibrating provides adequate results,asasalong assmall the used equationspistons. of acoustics are linear. To calculate sound media propagation as a set of small flat vibrating pistons. scattering properties reflecting boundaries, algorithm represents a geometric of Each wave elementoffalling onto such a pistonthe makes one radiate reflected sound in allmodel directiEach wavepropagation element falling such a piston makes one radiate sets reflected soundand in all directisound as construct aonto set of small flat vibrating pistons. ons andmedia makes it possible to an algorithm which accepts of sources reflecting ons and makes ityields possible to construct an algorithm whichpoints accepts sets of sources and reflecting surfaces. It also a field distribution over specified such that each source, Each wave element falling onto such a piston makes one radiate reflected sound in allprimary directisurfaces. It also yields a field distribution over specified points such that each source, primary or secondary, can be associated with an element of parallel execution and be managed via a list ons and makes it possible to construct an algorithm which accepts sets of sources and reflecting or secondary, can be associated with ana element of parallel execution and be managed via a list of polymorphic sources implementing task list. surfaces. It also yields a field distribution over specified points such that each source, primary of polymorphic sources implementing a task list. or secondary, be associated with an element ofofparallel execution and be via aused list The reportcan covers a mathematical formulation the problem, defines an managed object model The report the covers a mathematical formulation the problem, an object model used of polymorphic sources implementing a task some list. of to implement algorithm, and provides analysis of the defines algorithm in sequential and to implement algorithm, and provides some analysis of the algorithm in sequential and parallel forms. the The report covers a mathematical formulation of the problem, defines an object model used parallel forms. to2017 implement thePublished algorithm, and provides some analysis of the algorithm in sequential and © The Authors. by Elsevier B.V. Keywords: computer simulation, problem-oriented programming, high-performance computing, archiparallel forms. Keywords: under computer simulation, high-performance computing,Science archiPeer-review responsibility of the problem-oriented scientific committeeprogramming, of the International Conference on Computational tectural acoustics, parallel algorithms tectural acoustics, parallel algorithms Keywords: computer simulation, problem-oriented programming, high-performance computing, architectural acoustics, parallel algorithms
1 1
1877-0509 © 2017 The Authors. Published by Elsevier B.V. Peer-review under responsibility of the scientific committee of the International Conference on Computational Science 10.1016/j.procs.2017.05.161
1
1754
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
1
Introduction
Ray tracing is a commonly used approach to simulate wave fields in problem domains in which a wave nature of the field can be ignored, and, in particular, a given wavelength can be approximated to zero. This allows a fast rough simulation of the field which can be sufficient in a number of applications, including architectural acoustics. From the perspective of parallel computing, ray tracing allows independent analysis of each ray emitted by a source and thus can easily be implemented on a number of parallel devices, including graphic processors. The main disadvantage of the approach is the aforementioned neglect of wave properties during simulations, including diffractions and scatterings on media boundaries, although a set of directions in which the reflected wave has significant energy is rather wide. Therefore, the authors’ approach to use ray tracing and simultaneously take into account these effects is to consider reflections of sound in all directions, not just in one symmetrical to the incidence direction (with respect to a normal of the reflecting boundary). In order to provide such a mechanism the authors rely upon the Rayleigh’s integral over some flat surface dσ.
2
A mathematical model of the problem
Typical sound waves considered in architectural acoustics range from 2.1 centimeters to 6.8 meters in wavelengths (these wavelengths λ are obtained from a formula λ = c/f with sound velocity c ≈ 340m/s given for Normal Temperature and Pressure (NTP) and typical acoustical frequencies f from 50 Hz to 16 kHz). These wavelengths obviously cannot be neglected during simulations of the field in ordinary rooms of comparable sizes, so in the case of low and medium frequencies and comparable volumes of regions of field propagation media applications usually implement the Boundary Element Method (BEM) or the Finite Element Method (FEM). Implementations of these methods derived from some equilibrium equation of the field can yield rather precise results, especially the FEM. The downside is a heavy impact of a size of an element to the used computational resources (the element should be at most six times smaller than the considered wavelength [4]) and complications of handling absorbing boundary conditions often essential in the domain and which can involve heavy blocking synchronization of parallel threads. The presented below model has been developed and used by the authors in an attempt to mitigate these issues. The key element of the model is that it considers reflecting surface as a set of elementary vibrating pistons which reflect sound in multiply directions with a radiation pattern calculated via the Rayleigh’s integral: 1 ∂ϕ1 e−ikr2 Φ (M ) = dσ, (1) 2π ∂n r2 σ
∂ϕ1 ∂n
Here is a component of a vibrational speed of sound normal to the reflecting surface σ, r2 is a distance to some point M of interest, and k is a wave number. The integrand is a formula for an elementary point source. Therefore, the reflecting surface is represented as a set of such sources which cumulatively affect a state of the field in M . 1 When a flat surface is considered as a source by itself, a distribution of ∂ϕ ∂n is usually assumed constant throughout σ. But in the case of reflecting surface which vibrates as the result of an 2
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
interaction with a spherical wave, elements of the wave front approach to the surface with some delays in time which cause nonzero phase shifts in vibrations of each element dσ of the surface. With this nonuniformity taken into account, a primary point source radiates a spherical acoustic wave with its potential described by the formula (2), provided that an amplitude of the wave equals one unit. ei(ωt−kr1 ) , (2) ϕ1 = r1 where r1 is a distance from the primary source to the element dσ of the reflecting surface σ in (1), and ω is an angular frequency of sound. A differentiation of (2) with respect to the normal vector n yields: ∂ϕ1 ei(ωt−kr1 ) = ∇n ϕ1 = −ik cos γ ∂n r1
(3)
with γ being an angle between the vectors r1 (from the primary source to dσ) and n. Thus, provided that ν is a reflection coefficient dependent upon acoustical impedances of the media separated by the surface, reflection formula (1) takes the following form: i(ωt−k(r1 +r2 )) ikν e ϕ=− cos γdσ (4) 2π r1 r2 σ
2 −Z1 with ν = Z Z2 +Z1 , where Z1 and Z2 are the impedances with imaginary parts to be ignored due to a plain wave approximation for a sufficiently small dσ [5].
3
The object model
Below, there is a simplified domain object model presented here to define interfaces, in terms of which the algorithm is described below. Basically, the process of simulation is defined as a function of acoustical field distribution over given set of control points. An input parameter to the function is this set provided together with two sets of primary sources generating the field as specified by a user and the set of reflecting surfaces which constitute the medium. To implement such a function with an algorithm, the following top-level object model of the domain was built as follows (Fig. 1). The three input sets, together with some other spacial characteristics of the medium (e.g. dissipation), compose the encapsulating Medium class. The class aggregates all primary sources (instances of the PrimarySource class) of the medium and provides access to an enumerable set of them. Likewise, the class aggregates enumerable set of reflecting objects. The latter are implemented by the ReflectingObject class which similarly is represented as a set of flat media boundaries. They, in turn, can be discretized to sets elementary pistons dσ around points which appertain to the boundaries. Hence, the sets of these points are returned by the ReflectingElements method with respect to the considered frequency. Each of these reflecting elements dσ are associated with a secondary source implemented by the RayleighSource class. The class implements the integrand of (4) providing a polymorphic EmitTo method defined in the abstract Source base class. This gives an opportunity to treat all the sources of the field, both primary and secondary, uniformly. Also, the process is parametrized with a set of control points which accumulate characteristics of incoming sound and which are entities of the ControlPointSet enumerable class. The 3
1755
1756
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
Figure 1: Entities of the domain of acoustics
Figure 2: Visualization of sound distribution over points of a plane of the ControlPointSet class produced by the implementation below
AccumulateSound method of the class registers an effect of a specified source, on a given frequency, and at the region of a specified point (i.e. in the nearest point of the set). The approach used by the authors considers the set as a set of hyperplanes which are discretized onto a grid with the control points located in the nodes of the grid providing results such as one shown on 4
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
Fig. 2 with sound pressure distributed over such a plane.
4
The algorithm
The implementing algorithm indirectly exploits the formula (4) to calculate a field distribution over a given set of points. The field is generated by a set of primary sources which, in the turn, can produce sources associated with reflecting boundaries of the medium M . This requires some discretization process to yield boundary elements dσ for each of the reflecting surfaces constituting the medium. Therefore, each reflecting boundary is associated with a set of points appertaining to the boundary, so that each element dσ is constructed around each of the points. These sets are specified as a function of frequency, and each dσ is by itself considered a point source radiating reflected sound to the hemisphere of the incoming wave.
4.1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
A flow of the algorithm
Data: M :Medium Result: Distribution of sound on a set of plains in P S ← M .PrimarySources(); while |S | = 0 do forall si ∈ S do S ← S \ {si }; forall ω ∈ si .BandSet() do forall pP ∈ M .ControlPoints() do if M .RegistrationPredicate(s, pP ) then if ¬si .InShadowArea(M .ReflectingObjects(), pP ) then P .AccumulateSound(si , ωj , pP ) end end end forall r ∈ M .ReflectingObjects() do forall f ∈ r.ReflectingBoundaries() do forall pf ∈ f .ReflectingElements(ω) do if M .SecondarySourcePredicate(s, pf ) then if ¬si .InShadowArea(M .ReflectingObjects(), pf ) then s ← CreateRayleighSource(si , ω, pf ) ; S ← S ∪ {s } ; end end end end end end end end Algorithm 1: The Rayleigh integral based algorithm.
5
1757
1758
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
Figure 3: A parallel form of the algorithm
4.2
The parallel form
The presented algorithm can easily be parallelized as demonstrated by a graph of its parallel form (Fig. 3). The following possible fork-join pairs can exist. From the step 3 to 26, there is an outer loop which enumerates all the sources of the model, both primary and secondary, in a container of sources S which is synchronously modified during each iteration, i.e. when a new secondary source is obtained, or when the existing source is popped out of the container in order for the former to be processed. Thus, the set S of the above algorithm is essentially a task list with each one being specified by a polymorphic source object. Steps 5-25 constitute a branch of parallel execution for elements of the set W (si ) of bands. Provided that elements of the simulation for each ωj of W (si ) can be performed independently (which is the case), different iterations of the loop can be executed in parallel. Two inner loops 6-12 and 13-24 are mutually independent and therefore can be run in parallel as well, although this requires a support for MPMD-parallelism (Multiple Program Multiple Data). By itself, these loops can be parallelized over their iterations with the first one enumerating mutually independent control points which accumulate the simulation results, and the second one iterating over different reflecting objects of a read-only container of the medium M . The latter loop 13-24 performs a nested enumeration of points (i.e. elementary areas dσ in (4)) to construct secondary sources and insert them into the shared container S . Both of the nested loops can also be executed in parallel. This parallel form requires some blocking synchronization of threads of execution. Namely, on steps 4, 9, and 19, there are modifications of shared memory being made by the algorithm. These modifications should be guarded by a computationally heavy mutex synchronization. Also, on the step 2 a thread reads a state of the shared container S which should be protected by a condition variable synchronization with respect to the mutex synchronization on steps 4 and 19 [1], [2]. These points of access to shared data can result in significant performance degradation, especially in distributed computing environments, so divide and conquer based optimizations [3] might be necessary. In the implementation used by the authors task list S is distributed among nodes of a computer network with respect to their computational performance. The latter is measured as 6
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
a function of the number of logical CPU (Central Processor Unit) cores and a flag of the presence of graphical processor units capable for general-purpose computations (GPGPU). Each node works with its own instances of S , the set of control points, and the set of reflecting objects (still shared between local threads) with results in local control points combined with global ones after the local S is exhausted. The newly created tasks, i.e. secondary sources s on step 18, are posted to the global S directly. The distributed threads 3-26 run by the nodes fork at steps 5 and 15. If GPGPUs are present on a node, then, additionally, the loop 6-12 is executed on the GPGPUs. For more on the implementation and experimentally obtained effectiveness, see the sections 5 and 6 below.
4.3
The Complexity
In order to estimate a complexity of the algorithm, we associate an elementary step of execution with a sound wave element emitted by a source s ∈ S to some point p in the modeled space. A number of such elementary tasks created for a given source, not counting any derived sources produced during reflections of emitted sound, is a linear function
ns (si ) = |P | +
∀f ∈F
|f | · |W (si )|
(5)
which is determined by a number of control points, among which the sound distribution is sought, a numeric value of dσ in (4), and a number of bands in W (si ) = si .BandSet(). The latter is either determined by a discretization of the primary source frequency response or equals one for secondary sources defined by the reflected wave element. A more significant impact on the complexity is contributed by the total number of secondary sources produced during a simulation process. This number is defined by a value of I0 which is an imposed minimum of an obtained sound characteristic specified to limit the simulation on the step 15 of the algorithm. To obtain this number consider the average number nr of reflections which occur to a wave element emitted by a spherical source. An intensity of a sound wave scattered in space decreases as r−2 , and its pressure falls as r−1 . For simplicity consider a sound dissipation in the media to be neglectable which is often the case for small distances inside a room. Each reflection causes a ν-fold decrease of the obtained sound characteristic due to absorptions (ν is a reflection coefficient for a surface). Thus, before energy of an acoustical wave decreases to I0 , its elements are reflected nr times on average. Therefore, it can be given as nr = logν II0 , where I is the initial energy produced by the source si , and ν is an averaged reflection coefficient (6).
ν=
S (f ) · ν (f ) S (f )
∀f ∈F
(6)
∀f ∈F
where S(f ) is an area of a surface f , and ν(f ) is its reflection coefficient. Each source si produces N = |W (si )| · |f | secondary sources, then each of these secon∀f ∈F
dary sources produces the same number of tertiary sources and so on. Therefore, one primary n source produces the following sequence of secondary ones: 1 → N → N ·τ N → · · · → N ·(τ N ) r . Here τ is a constant correction needed to exclude a reflecting surface which a secondary source of consideration appertains to (e.g. 56 for box-shaped rooms). Thus, the sum of the members 7
1759
1760
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
of this geometrical progression is Ntotal = 1 +
n −1 N · 1 − (τ N ) r 1 − τN
.
(7)
Thus, there is an exponential dependency of a number of elementary simulations on a value of nr . Hence, the algorithm is feasible only for a small amount nr of reflections affecting a wave element emitted by a primary source. In an implementation used by the authors a value of nr is, at most, set to 4-5.
5
The implementation
The presented algorithm was implemented as a part of a distributed heterogeneous computer simulation system maintained by the authors. The system, which is designed to simulate wave physical fields, consists of the centralized structurally separated domain-independent subsystem (DIS) servicing the problem-oriented domain subsystem (DS). The former implements the service layer and is defined independently of a given problem domain, and the latter implements the domain layer. The DS is composed of implementations of simulations in specific problem domains (of wave field simulation) including one implementing the algorithm described in the report. The DIS provides general services to simulation implementations which include communication with a user and computing nodes of both subsystems, load balancing between nodes of simulation implementations, geometry modeling of the environment, and other ancillary services (Fig. 4). The geometry modeling is defined and managed separately by the DIS which conceptually considers the environment as a composition of a given medium, sources of field, reflecting boundaries, and a set of control points which accumulate the simulation results. The sources, the set of control points, and reflecting objects are all derived from a common class Object. The geometry model and its instances of the Object class are assigned their unique identifiers, associative arrays of domain-specific parameters (accessed using identifiers of simulation implementations), and string patterns specifying a set of identifiers of simulation implementations for which each element of the geometry model is defined. A chosen geometry model is sent to a specific simulation implementation along with corresponding domain parameters upon a request from a user to start a simulation process. When the request is received by a simulation implementation, the latter maps the geometry model (marshaled as either a binary package or an XML document) to the internal representation, such as one shown on Fig. 1. The load balancing is separated onto two levels. The first one is defined by the DIS and is used to distribute tasks among all nodes of the domain subsystem. Each task is specified with a set of its input parameters being stored by the load balancer with respect to the identifier of a simulation implementation from which the task has originated. Additionally, the load balancer implements priority queues of addresses of distributed nodes for each simulation implementation. The priority is specified by an integral number produced by the domain subsystem. The implementation of the algorithm defines this number as a function priority (node) = node.cpu cores () · (1 + (node.gpu cores () = 0 → 1) ∨ 0) .
(8)
Also, the implementation of acoustical fields simulation defines a similar load balancer, instances of which are deployed on every node which constitute the component. This load balancer is implemented as a queue of tasks and a pool of unoccupied threads (i.e. a total number of available logical CPU cores minus a number of threads performing the tasks) as 8
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
DIS -
Communication protocols Load balancing Geometry modeling Interface, data management, etc.
Simulation implementation 1 ... Simulation implementation n
Figure 4: The implementation on the top level
specified above. During the process of simulation, components of execution, which correspond to steps 3-16 of the algorithm and are specified with the aforementioned parameters of sources, are assigned to a free thread of the pool. If the pool is empty, the tasks are attempted to be forwarded to the load balancer of the DIS. If the latter is full (with respect to a specified maximum limit, in bytes, which herein is 128 GB) the parameters stay in the local queue. The implementation is deployed onto the following platform of five nodes. Four nodes use a two-processor configuration based on 12-core Intel Xeon X12 E5-2695v2 processors with HyperThreading and include CUDA Nvidia Tesla K40 GPU each. One of the four nodes is operated by Ubuntu Server 15.10. This node is used as a platform for the DIS. The other three are used for the deployment of the domain subsystem so that two of the clusters, which run on Windows Server 2012 R2, implement the presented algorithm. Additionally, the load balancer of the DIS uses a dedicated node to implement the data structures described above. This node has a RAID5 array of twelve hard drives and one Intel Xeon E5-2609v2 processor.
6
Analysis of the effectiveness of the algorithm
Functional effectiveness of the simulation is defined as a weighted sum of three basic indicators of effectiveness which are the time performance, space performance, and accuracy of the results. Time performance (table 1) is time measured on the user side from the moment of posting an asynchronous request to start a process of simulation for a given model (which has been posted to the DIS beforehand) until after the moment the user has been notified about the completion of the simulation. Space performance (table 1) includes memory consumptions by the request, its parameters, the geometry model, and maximum memory occupied by data structures of both load balancers described above. Accuracy (table 2) is determined for a given geometry model as an inversion of a standard deviation of intensity of the simulated acoustical field from values obtained from live measurements. To produce a numeric efficiency value, the results are also averaged for every control point and for every frequency and normalized (the result is not shown in the table). The experiments were performed in a hall of size 30 m×30 m×7.5 m with one source of sound emitting acoustical waves of frequencies 58 Hz to 18.5 kHz. The source has a cardioid radiation pattern, it was located in the position (14 m; 4.4 m; 1.2 m) and directed towards the farthest wall. Reflecting objects were represented in the geometry model by 164 polygonal surfaces. Measurements were obtained in points which appertain to a single plane parallel to the floor and 1.2 meters above it. The coefficient of attenuation in the medium of the hall was set to 1 (no attenuation) throughout the volume and spectrum. 9
1761
1762
Andrey A. Chusov et al. / Procedia Computer Science 108C (2017) 1753–1762 Parallel Acoustic Field Simulation Chusov, Statsenko, Lysenko, Cherkassova, Kuligin, Unru and Bernavskaya
Logical CPU cores involved Time performance, s Space consumptions, MiB
4 1.32E+6 10217
32 4.87E+5 10165
42 3.3E+5 10224
60 2.92E+5 10165
78 2.21E+5 9886
96 1.72E+5 9716
Table 1: Performance of the implementation Control point (2.8 m; 4.4 m; 1.2 m) (17.5 m; 4.4 m; 1.2 m) (2.8 m; 9 m; 1.2 m) (17.5 m; 9 m; 1.2 m) (2.8 m; 24.6 m; 1.2 m) (17.5 m; 24.6 m; 1.2 m)
Measured intensity, W m−2 6.92E-6 5.89E-4 8.32E-7 8.51E-5 6.17E-6 1.02E-5
Simulated intensity, W m−2 7.1E-6 5.95E-4 8.29E-7 7.54E-5 6.2E-6 9.77E-6
Ratio 1.03 1.01 1.0 0.89 1.01 0.88
Table 2: Accuracy on 250 Hz
7
Conclusion
The proposed algorithm of sound simulation takes into account sound scatterings during its reflections on surfaces constituting modeled medium of the acoustical field. This, on one hand, allows considering significant parts of secondary radiation back into space, but, at the same time, implies an exponential complexity over the average number of reflections which occur to a wave emitted by the field’s source. The algorithm can be easily parallelized on various numbers of parallel execution units, with most points of parallelization supporting the SPMD model. Experiments performed by the authors with the implementation of the algorithm demonstrate almost a full load of processors, parallel and distributed, used by the implementation and a hyperbolic dependency of a total time required by a single simulation upon a number of the processors. Nevertheless, this number affects the total simulation time only linearly and is significant only for a low average number of reflections which happen to different elements of an emitted acoustical wave.
References [1] James H. Anderson, Yong-Jik Kim, and Ted Herman. Shared-memory mutual exclusion: major research trends since 1986. Distributed Computing, 16(2):75–110, 2003. [2] Bj¨ orn B Brandenburg and James H Anderson. Reader-writer synchronization for shared-memory multiprocessor real-time systems. In Real-Time Systems, 2009. ECRTS’09. 21st Euromicro Conference on, pages 184–193, Los Alamitos, CA, USA, 2009. IEEE Computer Society. [3] Mattias Eriksson, Christoph Kessler, and Mikhail Chalabine. Load balancing of irregular parallel divide-and-conquer algorithms in group-SPMD programming environments. In Masters Thesis, PELAB, Linkpings Universitet, pages 313–322, 2006. [4] T. L odygowski and W. Sumelka. Limitations in application of finite element method in acoustic numerical simulation. Journal of Theoretical and Applied Mechanics, 44(4):849–865, 2006. [5] A. Moiola, R. Hiptmair, and I. Perugia. Plane wave approximation of homogeneous helmholtz solutions. Zeitschrift f¨ ur angewandte Mathematik und Physik, 62(5):809, 2011.
10
ScienceDirect
This space is reserved for the Procedia header, do not use it This space is reserved for the Procedia header, do not use it Procedia Computer Science 108C (2017) 1753–1762
This space is reserved for the Procedia header, do not use it
International Conference on Computational Science, ICCS 2017, 12-14 June 2017, Zurich, Switzerland
Parallel Parallel Acoustic Acoustic Field Field Simulation Simulation with with Respect Respect to to Scattering of Sound on Local Inhomogeneities Parallel Acoustic Field on Simulation with Respect to Scattering of Sound Local Inhomogeneities 1 1 of Sound on Local Inhomogeneities AndreyScattering A. Chusov11 , Lubov G. Statsenko 1 , Alexsey P. Lysenko1 , Sergey N.
Andrey1 ,A. Chusov , Lubov G.1 ,Statsenko , Alexsey P. Lysenko , Sergey N.2 1 Kuligin 1 Nina A. Cherkassova 1 Petr P. Unru 1 , and Maya V. Bernavskaya 2 1 1 1 Nina A. Cherkassova Petr P. Unru , and P. Maya V. Bernavskaya Kuligin Andrey ,A. Chusov , Lubov G. ,Statsenko , Alexsey Lysenko , Sergey N. 1 1 1 2 Far Eastern Federal University, Vladivostok, Russia Kuligin1 , Nina A. 1 Cherkassova , Petr P. Unru , and Maya V. Bernavskaya Far Eastern Federal University, Vladivostok, Russia 2 2 2
[email protected] [email protected] 1 Saint-Petersburg Polytechnic University, Saint-Petersburg, Far Eastern Federal University, Vladivostok, Russia Russia Saint-Petersburg Polytechnic University, Saint-Petersburg, Russia [email protected] [email protected] [email protected] Saint-Petersburg Polytechnic University, Saint-Petersburg, Russia [email protected]
Abstract Abstract The report presents a developed approach to simulation of acoustic fields in enclosed media. The report a developed to simulation of acoustic fields enclosed sources media. Abstract This methodpresents is based on the useapproach of Rayleigh’s integral for calculation of in secondary This method is based on the use of Rayleigh’s integral for calculation of secondary sources generated a wavea falling ontoapproach media boundaries. The is media. highly The reportby presents developed to simulation of implementing acoustic fields algorithm in enclosed generated by aimplies wave falling onto media boundaries. The implementing algorithm is highly parallelizable, loosely coupled parallel branches with only few points of inter-thread This method is based on the use of Rayleigh’s integral for calculation of secondary sources parallelizable, implies loosely coupled parallel branches with only few points of inter-thread communication. On the other hand, the algorithm isThe exponential upon an average generated by a wave falling onto media boundaries. implementing algorithm isnumber highly communication. On the other hand, the algorithm is exponential upon an average number of reflections which occur to a single wave element emitted by a primary source, although parallelizable, implies loosely coupled parallel branches with only few points of inter-thread of reflections which occurthis to a single can wave element emitted byprovide a primary source, although for practical applications reducedis enough to results with communication. On the othernumber hand, the be algorithm exponential uponaccurate an average number for practical applications this number can be reduced enough to provide accurate results with reasonable time and space of reflections which occur consumptions. to a single wave element emitted by a primary source, although reasonable time and space consumptions. for The practical applications thisis number reduced enough to provideofaccurate results proposed algorithm based oncan thebeapproximate superposition acoustical fieldswith and The proposed algorithm is based on the approximate superposition of acoustical and reasonable time and space consumptions. provides adequate results, as long as the used equations of acoustics are linear. To fields calculate provides adequate results, as long as the used equations of acoustics are linear. To calculate scattering properties of reflecting boundaries, the algorithm representsofa acoustical geometric fields modeland of The proposed algorithm is based on the approximate superposition scattering properties of reflecting boundaries, the algorithm represents a geometric model of sound media propagation set of flat vibrating provides adequate results,asasalong assmall the used equationspistons. of acoustics are linear. To calculate sound media propagation as a set of small flat vibrating pistons. scattering properties reflecting boundaries, algorithm represents a geometric of Each wave elementoffalling onto such a pistonthe makes one radiate reflected sound in allmodel directiEach wavepropagation element falling such a piston makes one radiate sets reflected soundand in all directisound as construct aonto set of small flat vibrating pistons. ons andmedia makes it possible to an algorithm which accepts of sources reflecting ons and makes ityields possible to construct an algorithm whichpoints accepts sets of sources and reflecting surfaces. It also a field distribution over specified such that each source, Each wave element falling onto such a piston makes one radiate reflected sound in allprimary directisurfaces. It also yields a field distribution over specified points such that each source, primary or secondary, can be associated with an element of parallel execution and be managed via a list ons and makes it possible to construct an algorithm which accepts sets of sources and reflecting or secondary, can be associated with ana element of parallel execution and be managed via a list of polymorphic sources implementing task list. surfaces. It also yields a field distribution over specified points such that each source, primary of polymorphic sources implementing a task list. or secondary, be associated with an element ofofparallel execution and be via aused list The reportcan covers a mathematical formulation the problem, defines an managed object model The report the covers a mathematical formulation the problem, an object model used of polymorphic sources implementing a task some list. of to implement algorithm, and provides analysis of the defines algorithm in sequential and to implement algorithm, and provides some analysis of the algorithm in sequential and parallel forms. the The report covers a mathematical formulation of the problem, defines an object model used parallel forms. to2017 implement thePublished algorithm, and provides some analysis of the algorithm in sequential and © The Authors. by Elsevier B.V. Keywords: computer simulation, problem-oriented programming, high-performance computing, archiparallel forms. Keywords: under computer simulation, high-performance computing,Science archiPeer-review responsibility of the problem-oriented scientific committeeprogramming, of the International Conference on Computational tectural acoustics, parallel algorithms tectural acoustics, parallel algorithms Keywords: computer simulation, problem-oriented programming, high-performance computing, architectural acoustics, parallel algorithms
1 1
1877-0509 © 2017 The Authors. Published by Elsevier B.V. Peer-review under responsibility of the scientific committee of the International Conference on Computational Science 10.1016/j.procs.2017.05.161
1
1754
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
1
Introduction
Ray tracing is a commonly used approach to simulate wave fields in problem domains in which a wave nature of the field can be ignored, and, in particular, a given wavelength can be approximated to zero. This allows a fast rough simulation of the field which can be sufficient in a number of applications, including architectural acoustics. From the perspective of parallel computing, ray tracing allows independent analysis of each ray emitted by a source and thus can easily be implemented on a number of parallel devices, including graphic processors. The main disadvantage of the approach is the aforementioned neglect of wave properties during simulations, including diffractions and scatterings on media boundaries, although a set of directions in which the reflected wave has significant energy is rather wide. Therefore, the authors’ approach to use ray tracing and simultaneously take into account these effects is to consider reflections of sound in all directions, not just in one symmetrical to the incidence direction (with respect to a normal of the reflecting boundary). In order to provide such a mechanism the authors rely upon the Rayleigh’s integral over some flat surface dσ.
2
A mathematical model of the problem
Typical sound waves considered in architectural acoustics range from 2.1 centimeters to 6.8 meters in wavelengths (these wavelengths λ are obtained from a formula λ = c/f with sound velocity c ≈ 340m/s given for Normal Temperature and Pressure (NTP) and typical acoustical frequencies f from 50 Hz to 16 kHz). These wavelengths obviously cannot be neglected during simulations of the field in ordinary rooms of comparable sizes, so in the case of low and medium frequencies and comparable volumes of regions of field propagation media applications usually implement the Boundary Element Method (BEM) or the Finite Element Method (FEM). Implementations of these methods derived from some equilibrium equation of the field can yield rather precise results, especially the FEM. The downside is a heavy impact of a size of an element to the used computational resources (the element should be at most six times smaller than the considered wavelength [4]) and complications of handling absorbing boundary conditions often essential in the domain and which can involve heavy blocking synchronization of parallel threads. The presented below model has been developed and used by the authors in an attempt to mitigate these issues. The key element of the model is that it considers reflecting surface as a set of elementary vibrating pistons which reflect sound in multiply directions with a radiation pattern calculated via the Rayleigh’s integral: 1 ∂ϕ1 e−ikr2 Φ (M ) = dσ, (1) 2π ∂n r2 σ
∂ϕ1 ∂n
Here is a component of a vibrational speed of sound normal to the reflecting surface σ, r2 is a distance to some point M of interest, and k is a wave number. The integrand is a formula for an elementary point source. Therefore, the reflecting surface is represented as a set of such sources which cumulatively affect a state of the field in M . 1 When a flat surface is considered as a source by itself, a distribution of ∂ϕ ∂n is usually assumed constant throughout σ. But in the case of reflecting surface which vibrates as the result of an 2
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
interaction with a spherical wave, elements of the wave front approach to the surface with some delays in time which cause nonzero phase shifts in vibrations of each element dσ of the surface. With this nonuniformity taken into account, a primary point source radiates a spherical acoustic wave with its potential described by the formula (2), provided that an amplitude of the wave equals one unit. ei(ωt−kr1 ) , (2) ϕ1 = r1 where r1 is a distance from the primary source to the element dσ of the reflecting surface σ in (1), and ω is an angular frequency of sound. A differentiation of (2) with respect to the normal vector n yields: ∂ϕ1 ei(ωt−kr1 ) = ∇n ϕ1 = −ik cos γ ∂n r1
(3)
with γ being an angle between the vectors r1 (from the primary source to dσ) and n. Thus, provided that ν is a reflection coefficient dependent upon acoustical impedances of the media separated by the surface, reflection formula (1) takes the following form: i(ωt−k(r1 +r2 )) ikν e ϕ=− cos γdσ (4) 2π r1 r2 σ
2 −Z1 with ν = Z Z2 +Z1 , where Z1 and Z2 are the impedances with imaginary parts to be ignored due to a plain wave approximation for a sufficiently small dσ [5].
3
The object model
Below, there is a simplified domain object model presented here to define interfaces, in terms of which the algorithm is described below. Basically, the process of simulation is defined as a function of acoustical field distribution over given set of control points. An input parameter to the function is this set provided together with two sets of primary sources generating the field as specified by a user and the set of reflecting surfaces which constitute the medium. To implement such a function with an algorithm, the following top-level object model of the domain was built as follows (Fig. 1). The three input sets, together with some other spacial characteristics of the medium (e.g. dissipation), compose the encapsulating Medium class. The class aggregates all primary sources (instances of the PrimarySource class) of the medium and provides access to an enumerable set of them. Likewise, the class aggregates enumerable set of reflecting objects. The latter are implemented by the ReflectingObject class which similarly is represented as a set of flat media boundaries. They, in turn, can be discretized to sets elementary pistons dσ around points which appertain to the boundaries. Hence, the sets of these points are returned by the ReflectingElements method with respect to the considered frequency. Each of these reflecting elements dσ are associated with a secondary source implemented by the RayleighSource class. The class implements the integrand of (4) providing a polymorphic EmitTo method defined in the abstract Source base class. This gives an opportunity to treat all the sources of the field, both primary and secondary, uniformly. Also, the process is parametrized with a set of control points which accumulate characteristics of incoming sound and which are entities of the ControlPointSet enumerable class. The 3
1755
1756
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
Figure 1: Entities of the domain of acoustics
Figure 2: Visualization of sound distribution over points of a plane of the ControlPointSet class produced by the implementation below
AccumulateSound method of the class registers an effect of a specified source, on a given frequency, and at the region of a specified point (i.e. in the nearest point of the set). The approach used by the authors considers the set as a set of hyperplanes which are discretized onto a grid with the control points located in the nodes of the grid providing results such as one shown on 4
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
Fig. 2 with sound pressure distributed over such a plane.
4
The algorithm
The implementing algorithm indirectly exploits the formula (4) to calculate a field distribution over a given set of points. The field is generated by a set of primary sources which, in the turn, can produce sources associated with reflecting boundaries of the medium M . This requires some discretization process to yield boundary elements dσ for each of the reflecting surfaces constituting the medium. Therefore, each reflecting boundary is associated with a set of points appertaining to the boundary, so that each element dσ is constructed around each of the points. These sets are specified as a function of frequency, and each dσ is by itself considered a point source radiating reflected sound to the hemisphere of the incoming wave.
4.1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
A flow of the algorithm
Data: M :Medium Result: Distribution of sound on a set of plains in P S ← M .PrimarySources(); while |S | = 0 do forall si ∈ S do S ← S \ {si }; forall ω ∈ si .BandSet() do forall pP ∈ M .ControlPoints() do if M .RegistrationPredicate(s, pP ) then if ¬si .InShadowArea(M .ReflectingObjects(), pP ) then P .AccumulateSound(si , ωj , pP ) end end end forall r ∈ M .ReflectingObjects() do forall f ∈ r.ReflectingBoundaries() do forall pf ∈ f .ReflectingElements(ω) do if M .SecondarySourcePredicate(s, pf ) then if ¬si .InShadowArea(M .ReflectingObjects(), pf ) then s ← CreateRayleighSource(si , ω, pf ) ; S ← S ∪ {s } ; end end end end end end end end Algorithm 1: The Rayleigh integral based algorithm.
5
1757
1758
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
Figure 3: A parallel form of the algorithm
4.2
The parallel form
The presented algorithm can easily be parallelized as demonstrated by a graph of its parallel form (Fig. 3). The following possible fork-join pairs can exist. From the step 3 to 26, there is an outer loop which enumerates all the sources of the model, both primary and secondary, in a container of sources S which is synchronously modified during each iteration, i.e. when a new secondary source is obtained, or when the existing source is popped out of the container in order for the former to be processed. Thus, the set S of the above algorithm is essentially a task list with each one being specified by a polymorphic source object. Steps 5-25 constitute a branch of parallel execution for elements of the set W (si ) of bands. Provided that elements of the simulation for each ωj of W (si ) can be performed independently (which is the case), different iterations of the loop can be executed in parallel. Two inner loops 6-12 and 13-24 are mutually independent and therefore can be run in parallel as well, although this requires a support for MPMD-parallelism (Multiple Program Multiple Data). By itself, these loops can be parallelized over their iterations with the first one enumerating mutually independent control points which accumulate the simulation results, and the second one iterating over different reflecting objects of a read-only container of the medium M . The latter loop 13-24 performs a nested enumeration of points (i.e. elementary areas dσ in (4)) to construct secondary sources and insert them into the shared container S . Both of the nested loops can also be executed in parallel. This parallel form requires some blocking synchronization of threads of execution. Namely, on steps 4, 9, and 19, there are modifications of shared memory being made by the algorithm. These modifications should be guarded by a computationally heavy mutex synchronization. Also, on the step 2 a thread reads a state of the shared container S which should be protected by a condition variable synchronization with respect to the mutex synchronization on steps 4 and 19 [1], [2]. These points of access to shared data can result in significant performance degradation, especially in distributed computing environments, so divide and conquer based optimizations [3] might be necessary. In the implementation used by the authors task list S is distributed among nodes of a computer network with respect to their computational performance. The latter is measured as 6
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
a function of the number of logical CPU (Central Processor Unit) cores and a flag of the presence of graphical processor units capable for general-purpose computations (GPGPU). Each node works with its own instances of S , the set of control points, and the set of reflecting objects (still shared between local threads) with results in local control points combined with global ones after the local S is exhausted. The newly created tasks, i.e. secondary sources s on step 18, are posted to the global S directly. The distributed threads 3-26 run by the nodes fork at steps 5 and 15. If GPGPUs are present on a node, then, additionally, the loop 6-12 is executed on the GPGPUs. For more on the implementation and experimentally obtained effectiveness, see the sections 5 and 6 below.
4.3
The Complexity
In order to estimate a complexity of the algorithm, we associate an elementary step of execution with a sound wave element emitted by a source s ∈ S to some point p in the modeled space. A number of such elementary tasks created for a given source, not counting any derived sources produced during reflections of emitted sound, is a linear function
ns (si ) = |P | +
∀f ∈F
|f | · |W (si )|
(5)
which is determined by a number of control points, among which the sound distribution is sought, a numeric value of dσ in (4), and a number of bands in W (si ) = si .BandSet(). The latter is either determined by a discretization of the primary source frequency response or equals one for secondary sources defined by the reflected wave element. A more significant impact on the complexity is contributed by the total number of secondary sources produced during a simulation process. This number is defined by a value of I0 which is an imposed minimum of an obtained sound characteristic specified to limit the simulation on the step 15 of the algorithm. To obtain this number consider the average number nr of reflections which occur to a wave element emitted by a spherical source. An intensity of a sound wave scattered in space decreases as r−2 , and its pressure falls as r−1 . For simplicity consider a sound dissipation in the media to be neglectable which is often the case for small distances inside a room. Each reflection causes a ν-fold decrease of the obtained sound characteristic due to absorptions (ν is a reflection coefficient for a surface). Thus, before energy of an acoustical wave decreases to I0 , its elements are reflected nr times on average. Therefore, it can be given as nr = logν II0 , where I is the initial energy produced by the source si , and ν is an averaged reflection coefficient (6).
ν=
S (f ) · ν (f ) S (f )
∀f ∈F
(6)
∀f ∈F
where S(f ) is an area of a surface f , and ν(f ) is its reflection coefficient. Each source si produces N = |W (si )| · |f | secondary sources, then each of these secon∀f ∈F
dary sources produces the same number of tertiary sources and so on. Therefore, one primary n source produces the following sequence of secondary ones: 1 → N → N ·τ N → · · · → N ·(τ N ) r . Here τ is a constant correction needed to exclude a reflecting surface which a secondary source of consideration appertains to (e.g. 56 for box-shaped rooms). Thus, the sum of the members 7
1759
1760
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
of this geometrical progression is Ntotal = 1 +
n −1 N · 1 − (τ N ) r 1 − τN
.
(7)
Thus, there is an exponential dependency of a number of elementary simulations on a value of nr . Hence, the algorithm is feasible only for a small amount nr of reflections affecting a wave element emitted by a primary source. In an implementation used by the authors a value of nr is, at most, set to 4-5.
5
The implementation
The presented algorithm was implemented as a part of a distributed heterogeneous computer simulation system maintained by the authors. The system, which is designed to simulate wave physical fields, consists of the centralized structurally separated domain-independent subsystem (DIS) servicing the problem-oriented domain subsystem (DS). The former implements the service layer and is defined independently of a given problem domain, and the latter implements the domain layer. The DS is composed of implementations of simulations in specific problem domains (of wave field simulation) including one implementing the algorithm described in the report. The DIS provides general services to simulation implementations which include communication with a user and computing nodes of both subsystems, load balancing between nodes of simulation implementations, geometry modeling of the environment, and other ancillary services (Fig. 4). The geometry modeling is defined and managed separately by the DIS which conceptually considers the environment as a composition of a given medium, sources of field, reflecting boundaries, and a set of control points which accumulate the simulation results. The sources, the set of control points, and reflecting objects are all derived from a common class Object. The geometry model and its instances of the Object class are assigned their unique identifiers, associative arrays of domain-specific parameters (accessed using identifiers of simulation implementations), and string patterns specifying a set of identifiers of simulation implementations for which each element of the geometry model is defined. A chosen geometry model is sent to a specific simulation implementation along with corresponding domain parameters upon a request from a user to start a simulation process. When the request is received by a simulation implementation, the latter maps the geometry model (marshaled as either a binary package or an XML document) to the internal representation, such as one shown on Fig. 1. The load balancing is separated onto two levels. The first one is defined by the DIS and is used to distribute tasks among all nodes of the domain subsystem. Each task is specified with a set of its input parameters being stored by the load balancer with respect to the identifier of a simulation implementation from which the task has originated. Additionally, the load balancer implements priority queues of addresses of distributed nodes for each simulation implementation. The priority is specified by an integral number produced by the domain subsystem. The implementation of the algorithm defines this number as a function priority (node) = node.cpu cores () · (1 + (node.gpu cores () = 0 → 1) ∨ 0) .
(8)
Also, the implementation of acoustical fields simulation defines a similar load balancer, instances of which are deployed on every node which constitute the component. This load balancer is implemented as a queue of tasks and a pool of unoccupied threads (i.e. a total number of available logical CPU cores minus a number of threads performing the tasks) as 8
Parallel Acoustic Field Simulation Lysenko, Cherkassova, Unru and Bernavskaya Andrey A.Chusov, Chusov etStatsenko, al. / Procedia Computer Science 108CKuligin, (2017) 1753–1762
DIS -
Communication protocols Load balancing Geometry modeling Interface, data management, etc.
Simulation implementation 1 ... Simulation implementation n
Figure 4: The implementation on the top level
specified above. During the process of simulation, components of execution, which correspond to steps 3-16 of the algorithm and are specified with the aforementioned parameters of sources, are assigned to a free thread of the pool. If the pool is empty, the tasks are attempted to be forwarded to the load balancer of the DIS. If the latter is full (with respect to a specified maximum limit, in bytes, which herein is 128 GB) the parameters stay in the local queue. The implementation is deployed onto the following platform of five nodes. Four nodes use a two-processor configuration based on 12-core Intel Xeon X12 E5-2695v2 processors with HyperThreading and include CUDA Nvidia Tesla K40 GPU each. One of the four nodes is operated by Ubuntu Server 15.10. This node is used as a platform for the DIS. The other three are used for the deployment of the domain subsystem so that two of the clusters, which run on Windows Server 2012 R2, implement the presented algorithm. Additionally, the load balancer of the DIS uses a dedicated node to implement the data structures described above. This node has a RAID5 array of twelve hard drives and one Intel Xeon E5-2609v2 processor.
6
Analysis of the effectiveness of the algorithm
Functional effectiveness of the simulation is defined as a weighted sum of three basic indicators of effectiveness which are the time performance, space performance, and accuracy of the results. Time performance (table 1) is time measured on the user side from the moment of posting an asynchronous request to start a process of simulation for a given model (which has been posted to the DIS beforehand) until after the moment the user has been notified about the completion of the simulation. Space performance (table 1) includes memory consumptions by the request, its parameters, the geometry model, and maximum memory occupied by data structures of both load balancers described above. Accuracy (table 2) is determined for a given geometry model as an inversion of a standard deviation of intensity of the simulated acoustical field from values obtained from live measurements. To produce a numeric efficiency value, the results are also averaged for every control point and for every frequency and normalized (the result is not shown in the table). The experiments were performed in a hall of size 30 m×30 m×7.5 m with one source of sound emitting acoustical waves of frequencies 58 Hz to 18.5 kHz. The source has a cardioid radiation pattern, it was located in the position (14 m; 4.4 m; 1.2 m) and directed towards the farthest wall. Reflecting objects were represented in the geometry model by 164 polygonal surfaces. Measurements were obtained in points which appertain to a single plane parallel to the floor and 1.2 meters above it. The coefficient of attenuation in the medium of the hall was set to 1 (no attenuation) throughout the volume and spectrum. 9
1761
1762
Andrey A. Chusov et al. / Procedia Computer Science 108C (2017) 1753–1762 Parallel Acoustic Field Simulation Chusov, Statsenko, Lysenko, Cherkassova, Kuligin, Unru and Bernavskaya
Logical CPU cores involved Time performance, s Space consumptions, MiB
4 1.32E+6 10217
32 4.87E+5 10165
42 3.3E+5 10224
60 2.92E+5 10165
78 2.21E+5 9886
96 1.72E+5 9716
Table 1: Performance of the implementation Control point (2.8 m; 4.4 m; 1.2 m) (17.5 m; 4.4 m; 1.2 m) (2.8 m; 9 m; 1.2 m) (17.5 m; 9 m; 1.2 m) (2.8 m; 24.6 m; 1.2 m) (17.5 m; 24.6 m; 1.2 m)
Measured intensity, W m−2 6.92E-6 5.89E-4 8.32E-7 8.51E-5 6.17E-6 1.02E-5
Simulated intensity, W m−2 7.1E-6 5.95E-4 8.29E-7 7.54E-5 6.2E-6 9.77E-6
Ratio 1.03 1.01 1.0 0.89 1.01 0.88
Table 2: Accuracy on 250 Hz
7
Conclusion
The proposed algorithm of sound simulation takes into account sound scatterings during its reflections on surfaces constituting modeled medium of the acoustical field. This, on one hand, allows considering significant parts of secondary radiation back into space, but, at the same time, implies an exponential complexity over the average number of reflections which occur to a wave emitted by the field’s source. The algorithm can be easily parallelized on various numbers of parallel execution units, with most points of parallelization supporting the SPMD model. Experiments performed by the authors with the implementation of the algorithm demonstrate almost a full load of processors, parallel and distributed, used by the implementation and a hyperbolic dependency of a total time required by a single simulation upon a number of the processors. Nevertheless, this number affects the total simulation time only linearly and is significant only for a low average number of reflections which happen to different elements of an emitted acoustical wave.
References [1] James H. Anderson, Yong-Jik Kim, and Ted Herman. Shared-memory mutual exclusion: major research trends since 1986. Distributed Computing, 16(2):75–110, 2003. [2] Bj¨ orn B Brandenburg and James H Anderson. Reader-writer synchronization for shared-memory multiprocessor real-time systems. In Real-Time Systems, 2009. ECRTS’09. 21st Euromicro Conference on, pages 184–193, Los Alamitos, CA, USA, 2009. IEEE Computer Society. [3] Mattias Eriksson, Christoph Kessler, and Mikhail Chalabine. Load balancing of irregular parallel divide-and-conquer algorithms in group-SPMD programming environments. In Masters Thesis, PELAB, Linkpings Universitet, pages 313–322, 2006. [4] T. L odygowski and W. Sumelka. Limitations in application of finite element method in acoustic numerical simulation. Journal of Theoretical and Applied Mechanics, 44(4):849–865, 2006. [5] A. Moiola, R. Hiptmair, and I. Perugia. Plane wave approximation of homogeneous helmholtz solutions. Zeitschrift f¨ ur angewandte Mathematik und Physik, 62(5):809, 2011.
10