Realizations of the generalized adaptive cross approximation in an acoustic time domain boundary element method

In acoustics, the boundary element method (BEM) is much more common compared to elasticity. This is driven by the applications, which are in acoustics very often radiation problems causing trouble in finite element or finite volume methods. Nevertheless, for a time domain calculation, still an efficient BE formulations is lacking. We consider the time domain BEM for the homogeneous wave equation with vanishing initial conditions and given Dirichlet boundary conditions. The generalized convolution quadrature method (gCQ) developed by Lopez‐Fernandez and Sauter is used for the temporal discretisation. The spatial discretisation is done classically. Essentially, the gCQ requires to establish boundary element matrices of the corresponding elliptic problem in Laplace domain at several complex frequencies. Consequently, an array of system matrices is obtained. This array of system matrices can be interpreted as a three‐dimensional array of data, which should be approximated by a data‐sparse representation. The adaptive cross approximation (ACA) can be generalized to handle these three‐dimensional data arrays. Adaptively, the rank of the three‐dimensional data array is increased until a prescribed accuracy is obtained. On a pure algebraic level, it is decided whether a low‐rank approximation of the three‐dimensional data array is close enough to the original matrix. Hierarchical matrices are used in the two spatial dimensions and the third dimension are the complex frequencies. Hence, the algorithm makes not only a data sparse approximation in the two spatial dimensions but detects adaptively how much frequencies are necessary for which matrix block.


INTRODUCTION
Time domain boundary element methods (BEM) are in acoustics an alternative method with special advantages in case of scattering.The basis are boundary integral equations as counterpart to the governing hyperbolic partial differential equation.The mathematical theory goes back to the beginning of the last century and can be found in the textbook by Sayas [1].The first realisation in time domain is originated by Mansur [2] in the 80th of the last century.Beside this approach, utilising a space-time ansatz the convolution quadrature-based formulations exist (see, e.g.Schanz [3]).The convolution quadrature method proposed by Lubich [4] allows to establish more stable methods.The generalization of this seminal technique to variable time step sizes has been proposed by López-Fernández and Sauter [5,6] and is called generalized convolution quadrature method (gCQ).Despite several applications of the convolution quadrature-based time domain BEM (see, e.g.Refs.[7,8]), an efficient formulation is missing.All BE formulations either for elliptic as well as for hyperbolic problems need so-called fast methods to handle real-world problems.If the gCQ is used, essentially, a three-dimensional data array has to be efficiently computed and stored.Here, the generalized Adaptive Cross Approximation (3D-ACA) will be used to obtain a data sparse formulation.This technique is a generalization of ACA [9] and is proposed by Bebendorf [10,11].It is based on a Tucker decomposition [12] and can be traced back to the group of Tyrtyshnikov [13].The application of the 3D-ACA can as well be found in Seibel [14].Different to the approach here, it has been applied on a standard convolution quadrature-based formulation and  2 -matrices are used.Here, it is applied on the indirect BE formulation using the single layer potential and the gCQ as time discretisation.For the spatial compression, usual -matrices are used with ACA.

PROBLEM SETTING
Let Ω ⊂ ℝ 3 be a bounded Lipschitz domain and Γ ∶= Ω its boundary with the outward normal .Denoting the wave speed of the fluid with , the acoustic wave propagation is governed by with the end time  > 0. The Dirichlet trace is denoted with  0 and its application gives the boundary data   .It is well known that the single layer potential is a solution of Equation (1).The unknown density function  is determined by the boundary integral equation The fundamental solution corresponding to the PDE (1a) is denoted by a capital letter ( − ,  − ).

Spatial discretisation:
The boundary Γ is discretised with elements resulting in an approximation Γ ℎ = ⋃  =1   , which is the union of geometrical elements.  denote boundary elements, here linear surface triangles, and their total number is .Next, finite element bases on the boundary Γ ℎ are used to construct the approximation space  = Span{ 1 ,  2 , … ,   }.The unknown density  is approximated by a linear combination of functions in Note, the coefficients   () are still functions of the continuous time .In the following, the shape functions   will be chosen constant discontinuous, that is,  =  holds.This is in accordance with the necessary function spaces for the integral equation, which are not discussed here (see Sayas [1]).Inserting these spatial shape functions in the boundary integral equation ( 3) and applying a collocation method results in the semi-discrete equation system  *  ℎ =   . ( The given boundary data are approximated by the same shape functions.Note, in Equation ( 5) in the vectors, the nodal values are collected and the matrix  contains the respective values of the spatially integrated fundamental solution.

Temporal discretisation:
Discretising the time in  not necessarily constant time steps Δ  , that is, [0, ] = [0,  1 ,  2 , … ,   ], Δ  =   −  −1 ,  = 1, 2, … , , the gCQ method [5,6] can be applied.Let us assume an A-and L-stable Runge-Kutta method given by its Butcher tableau     with  ∈ ℝ × , ,  ∈ ℝ  and  is the number of stages.The stability assumptions require that    −1 = (0, 0, … , 1) holds.With the vector 1 = (1, 1, … , 1) T of size , the semi-discrete boundary integral equation ( 5) can be approximated by The notation ()  indicates the discrete value of the respective function/vector at   and V indicates that the Laplace transform of the integral kernel is used in the single layer operator.Details on the algorithm, especially the used parameters, can be found in Lopez-Fernandez and Sauter [15].The brief summary of the solution procedure is given next: with implicit assumption of zero initial condition.• For all steps  = 2, … , , the algorithm has two steps 1. Update the solution vector  −1 at all integration points for  = 1, … ,   with the number of integration points   .2. Compute the solution of the integral equation at the actual time step Note the dependency of  −1 , which is the influence of the time history.
Essentially, this algorithm requires the evaluation of the integral kernel at   points   , which are complex frequencies.Consequently, we get an array of system matrices  ∈ ℂ ××  .This array of system matrices can be interpreted as a three-dimensional array of data which will be approximated by a data-sparse representation based on 3D-ACA.
A L G O R I T H M 1 Pseudo code of 3D-ACA (taken from Seibel [14]).

GENERALIZED ADAPTIVE CROSS APPROXIMATION: 3D-ACA
An approximation of a three-dimensional array of data or a tensor of third order  ∈ ℂ ××  in terms of a low-rank approximation has been proposed in Bebendorf et al. [11] and is referred to as a generalization of adaptive cross approximation or also denoted 3D-ACA.In this approach, the 3D array of data to be approximated is generally generated by defining the outer product with  ∈ ℂ × ,  ∈ ℂ   .The matrix  corresponds to the spatial discretization of the single layer potential at a specific frequency   , that is, its elements are determined at all collocation points   and for all elements .This matrix will be called face or slice and may be computed as a dense matrix or with a suitable hierarchical partition of the mesh as matrix.Within this matrix, different low-rank approximations can be applied, where here ACA [9] is used.The vector  collects selected elements of  at a set of frequencies determined by the gCQ.The latter would amount in   entries and the same amount of slices.The 3D-ACA on the one hand approximates, as above mentioned, the slices with low rank matrices and, more importantly, the used frequencies adaptively until a prescribed precision  compared to  is obtained.The latter is measured with the Frobenius norm.The basic concept of this approach is sketched in Algorithm 1.
The stopping criterion requires a norm evaluation of  () .As stated in Seibel [14], under the assumption of monotonicity, the norm can be computed iteratively by ) Like in the 'normal' ACA, the crosses, which are here three-dimensional crosses, can be selected arbitrarily.However, selecting the maximum element ensures convergence.For details, see Bebendorf et al. [11] or the application using  2matrices in the slices [14].But in the current case, there is an open point regarding the determination of a maximal element in an -matrix with low rank blocks.Certainly, it is no option to compute the matrix entry of face  by multiplying the low-rank matrices.Instead, the index  or  is defined by the maximum of the absolute row sum of the low rank matrices, respectively.To get the idea, let us assume a block  within the face approximated by the truncated singular value decomposition  ≈   = Σ  .In this situation, the low-rank matrices  and  consist of orthonormal columns and, hence, the index  of the pivot position can be determined with The analogous way is used also for index .This concept is based on the fact that the columns of  and  are orthonormal.
To adopt this to an ACA approximated block, a QR-decomposition of the low-rank matrix  and  is performed and, second, the SVD is applied on the smaller inner matrix      .Hence, the same structure and properties as for an SVD approximated block is obtained and Equation ( 11) can be used to determine the pivot index.The sketched technique is known as recompression in the ACA literature, which as a side effect improves the compression rates.

NUMERICAL STUDIES
The presented 3D-ACA is used in Equation ( 6) to obtain a data-sparse realisation of the demanding time domain BE formulation.There are two possibilities how to apply the 3D-ACA.Either it can be applied on each matrix block, which requires that for the different frequencies, the same hierarchical structure is used.This holds here as a pure geometrical admissibility criterium for the -matrix has been selected.The other choice is to apply 3D-ACA on the -matrix as a whole.Both versions will be considered in the following tests.
The proposed methodology is tested by computing the response of a unit cube [−0.5, 0.5] 3 due to an excitation at  = (0.8, 0.2, 0.3)  by a smooth impulse As mentioned above, linear triangles are used and constant shape functions for the unknown density.The underlying time stepping method has been the 2-stage Radau IIA.The approximation within the faces has been chosen to be   = 10 −4 … 10 −8 for the different refinement levels of the spatial discretisation.The latter starts with a meshsize ℎ = 0.5 m.In each refinement level, the mesh size is halved as well as the time step size.The ratio of time steps to mesh size was kept constant with 0.7.The precision of the method with respect to the frequencies, that is, the  in Algorithm 1 was selected as  = 100 *   .
For the first comparison, the error in space and time is observed at a set of inner points in Figure 1.With the selected constant shape functions, the spatial discretisation is the limiting factor.The 2-stage Radau IIA would allow higher orders.Hence, the quadratic convergence is observed for both, the dense solution as well as for the solution using 3D-ACA.Three different data sparse solutions are presented: having in the faces dense matrices denoted by 'none', having -matrices with ACA denoted 'ACA', and the same but with the recompression algorithm not only used for the determination of the pivot element as discussed above but also for the stored matrix blocks.The later is denoted 'RE-ACA'.The presented results in Figure 1 are those with 3D-ACA in the blocks.Certainly, the same results are obtained if 3D-ACA is applied on the whole matrix.After this check of a correct implementation and showing that the compression does not spoil the results, the efficiency with respect to the used frequencies, that is, the rank in the third matrix dimension is studied.In Figure 2, the number  of used frequencies is compared to a dense calculation.Obviously, very few frequencies are necessary to obtain the same error.The application of 3D-ACA on the matrix blocks seem to be preferred as even the mean value of the different matrix blocks is significantly smaller then applying the 3D-ACA on the whole matrix.This effect is mainly due to the stopping criterium using the norm (10) over the whole matrix and not only within one block.If the faces are also compressed as displayed in Figures 3 and 4, the tendency is the same.However, the maximum number of frequencies for different blocks become much more volatile.The latter is caused by the stopping criterion.Essentially, a continuous decrease of the approximation error due to 3D-ACA is assumed but the numerical realization shows oscillations at small error levels.This causes that the stopping criterion is sometimes nearly reached but not good enough.Overall, it can be concluded that the application of 3D-ACA on the blocks in the -matrix is the preferable version.
The above statement is confirmed by the obtained compression rates.These are the used memory compared to a dense computation of the overall algorithm.These numbers are plotted in Figure 5 versus the refinement.
The tendency of the block version is much better than that of the whole matrix.It must be remarked that the main compression steams from the reduction of frequencies and not from any compression within the faces.Hence, one or two frequencies more result in dramatically worse compression rates.But overall, the algorithm is very efficient because the storage is reduced up to 0.6%.

CONCLUSION
Time domain boundary element formulations have due to the convolution in time a very high demand on storage.For realworld problems, data-sparse methods are necessary.For the indirect BE formulation of a Dirichlet problem in acoustics, a data-sparse method is presented.It utilises the standard spatial discretization and the gCQ in time.To reduce the storage of the necessary three-dimensional data array, the generalized adaptive cross approximation (3D-ACA) is applied.The proposed method shows that a significant storage reduction is possible, up to 0.6% in the selected numerical example.Further, it is shown that 3D-ACA should be applied on each block of an hierarchical matrix representation of the system matrix.This requires the same hierarchical partitioning for all frequencies but is much more efficient than applying the 3D-ACA on the whole matrix.

F I G U R E 1
Convergence rate at inner points: dense and 3D-ACA results.F I G U R E 2 Necessary frequencies: no compression in the face.F I G U R E 3 Necessary frequencies: faces with ACA.F I G U R E 4 Necessary frequencies: faces with recompressed ACA.F I G U R E 5 Compression of the data-sparse representation.
Financial support by the joint DFG/FWF Collaborative Research Centre CREATOR (CRC -TRR361/F90) is gratefully acknowledged.O R C I D Martin Schanz https://orcid.org/0000-0002-6177-8751R E F E R E N C E S