Time-domain
Yes, this is time-domain, not frequency-domain!
Fast
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:
T. Takahashi, M. Tanigawa, N. Miyazawa, "An enhancement of the fast time-domain boundary element method for the three-dimensional wave equation", Computer Physics Communications, Vol. 271, 108229, 2022. (arXiv)
Memory-efficient
The FMM can result in a less amount of memory than the conventional TDBEM. This feature is useful when you solve large-scale problems.
Parallelizable /Vectorizable
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.
General-purpose
The aXustica can handle arbitrary-shaped surface(s) by virtue of triangle mesh. In the case of exterior problems, where the computational domain spreads infinitely, we do NOT have to consider any absorbing boundary condition, e.g. the PML for the case of the FDTD method. 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.
Stable
This is a result of the recent progress on the time stabilization of the TDBEM; see the paper in arXiv. Specifically, the Burton-Miller-type boundary integral equation is employed.
Acoustic analysis with a Nagoya university model
The buildings EI, IB, and No.2 of Engineering school of Nagoya university are considered as scatteres. The sound source is shown as a small sphere. The ground is considered by considering the mirror image. All the surfaces of the buildings (as well as the ground) are sound-hard. The map on the RHS is from the web site of the university. (2025/2/19)
One million scale simulation for a terrain field
The aXustica was applied to analyze the sound propagation over a terrain field. It should be noted that the field is represented by an open surface, i.e., the field is a square region of 0<x<L and 0<x<L, where L is 6800 m. Accordingly, since a point source, which is a sinusoidal wave with the wavelength of about 340 m (corresponding to about 1 Hz when the wave velocity c is 340 m/s), is given at x=L, y=L/2, and z=L/10, the simulation is valid until the sound arrives at the border of the open surface. Then, the final time T is determined as (L/2)/c = 20 s. The model consists of 1,045,458 triangles and available from TURBOSQUID. The time-step size Dt was given as 0.05 s and thus the number of time steps was determined as T/Dt=200. The analysis was done with two workstations. The total computation time was 14,981 s. Each workstation spent about 320 GB memory. (2024/12/25)
Top view of the model.
Mesh around the highest mountain that locates near the center of the surface.
Animation for viewing the model.
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.
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.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. 1 Initial shape
Fig. 2 Optimal (final) shape
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:
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.
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.
The 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 figures) 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.
Surface mesh for the dinosaur model.
Distribution of the sound pressure on the model.
Comparison of the fast TDBEM ("FMM") with the conventional TDBEM ("Conv."). The numbers 6, 8, 10, and 12 for "FMM" denotes the precision parameters; larger is more accurate but slower.