aXustica

Time-domain boundary element method  

for the 3D wave equation

Yes, this is time-domain, not frequency-domain!

The method is accelerated by a variant of the fast multipole method (FMM), which is similar to the planewave time-domain (PWTD) algorithm but more algebraic rather than analytic because it exploits the interpolation-based FMM. The computational complexity of the  fast time-domain boundary element method (TDBEM) is estimated as O(Ns^(1+f)*Nt), while that of the conventional TDBEM is O(Ns^2*Nt), where the number f is 1/3 or 1/2 (according to the spatial distribution of the boundary elements) and Ns and Nt denote the spatial and temporal DOFs, respectively. The details are described in the papers:

The FMM can result in a less amount of memory than the conventional TDBEM. This feature is useful when you solve large-scale problems.

The aXustica can be parallelized with OpenMP on a shared-memory machine. Thus, you can fully utilize all the computing cores on your PC. Moreover, time-consuming parts can be vectorized with AVX512 instructions. Further, the aXustica can be parallelized on a distributed computing system with the help of the OpenMPI.

The aXustica can handle arbitrary-shaped surface(s) in a computational domain. In the case of external problems, where the computational domain spreads infinitely, we do NOT have to consider any absorbing boundary condition (ABC; e.g. the PML in FDTD method) on the boundary of the domain. You can give Dirichlet and/or Neumann boundary condition on such one or more closed surfaces. In addition, you can program the incident field that you want, e.g., planewave, point source. 

This is a result of the recent progress on the time stabilization of the TDBEM; see the paper in arXiv. The stabilization is indispensable when one solves external Dirichlet problems, while the stability is not serious for external Neumann problems, which have been mainly considered so far. 

Performance: FMM vs Conventional method

The following table shows a numerical result of an exterior Neumann problem, where the scatter is discretized with 13,240 boundary elements (which models a dinosaur in the RHS figure) and the number of time steps is 480. In the table, "Conv." stands for the conventional TDBEM and the number following "FMM" denotes the precision parameter; larger is more accurate . First of all, the FMM is >10x faster than the conventional method, where both methods are parallelized with 20 cores on dual Xeon E5-2687W. In addition, you can see that the precision parameter can control the trade-off between speed (as well as memory) and accuracy.

Verification: Scattering by a sphere

It is possible to compute the sound pressure at any point. This computation is also accelerated by the FMM. This example demonstrates such a calculation. Specifically, we considered a scattering problem due to a acoustically-soft sphere when a continuous pulse is impinged to it. The sound pressure on a cross section (width 4 and height 3) including the center of the sphere is available as the animation. The number of boundary elements and interior points are 5120 and 4428, respectively, while the number of time steps is 400. The computation time was about 292 second with the above-mentioned 20 cores. 

Sound pressure at the final (399) time step. This picture was created by Gnuplot in a stupid way. So the jaggy boundary profile is not the actual one, which consists of 5120 triangles. It should be noted that the computational domain is INFINITE and the sound pressure can be computed at any point; the present FMM-accerated BEM inherits all the advantages of the BEM!

Simulation: Room Acousitics

As another example, the animation of Figure 18 in the JCP paper is found here. You can see the sound generated in the RHS room at t=0 propagates to the LHS room and radiates outside through the windows and door. In the animation, the outer surfaces are removed to see the inside of the house.

Shape optimization

We exploited the aXustica to solve a shape optimization problem. Specifically, we determined the shapes of two rigid cubes so that the temporal-integral of the squared sound pressure at a predefined observation point is maximized; See Fig. 1. The incident wave, which is a planewave pulse, is impinged from the front (i.e. -x side) in Fig. 1. 

Fig. 1 Initial shape

Fig. 2 Optimal (final)  shape

Fig.2 shows the optimal shape, by which the value of the present objective functional J is increased from 0.14 to 0.63; See Fig. 3. In addition, Fig. 4 shows the profile of the sound pressure p at the observation point.

Fig.3 History of ojbective functional J.

Fig. 4 Profiles of the sound pressure at the observation point.

The following animations show the distribution of the sound pressure p on the plane (z=0.5),  which includes the observation point. 

Initial shape

Optimal shape

The present analysis was performed with the former graduate student N. Miyazawa at Nagoya University.

The details of the shape optimization with aXustica is found in the paper:

T. Takahashi, N. Miyazawa, M. Tanigawa, "A three-dimensional shape optimization for transient acoustic scattering problems using the time-domain boundary element method", International Journal for Numerical Methods in Engineering, Vol. 124, pp.482-512, 2023.

Hybrid-parallelization on a memory-distributed system

An OpenMP/OpenMPI hybrid parallelization is implemented to run the aXustica on a memory-distributed computing system. The following figures show the result on a different number of boundary elements with one to four homogeneous computing nodes. The code is vectorized with the AVX512 instructions. (2023/02/20)

Total computation time with 1, 2, 3, and 4 computing nodes.

NOTE: The number of time steps is not a constant but increases as the number of boundary elements increases.

Speedup of using 2, 3, and 4 nodes against using 1 node.