Faisal Amlani

Stéphanie Chaillat

Adrien Loseille

An efficient preconditioner for adaptive Fast Multipole accelerated Boundary Element Methods to model time-harmonic 3D wave propagation

Keywords: Boundary Element Method, Fast Multipole Method, Anisotropic meshes, Preconditioning, Hierarchical matrices

This paper presents an efficient algebraic preconditioner to speed up the convergence of Fast Multipole accelerated Boundary Element Methods (FM-BEMs) in the context of time-harmonic 3D wave propagation problems and in particular the case of highly non-uniform discretizations. Such configurations are produced by a recently-developed anisotropic mesh adaptation procedure that is independent of partial differential equation and integral equation. The new preconditioning methodology exploits a complement between fast BEMs by using two nested GMRES algorithms and rapid matrix-vector calculations. The fast inner iterations are evaluated by a coarse hierarchical matrix (H-matrix) representation of the BEM system. These inner iterations produce a preconditioner for FM-BEM solvers. It drastically reduces the number of outer GMRES iterations. Numerical experiments demonstrate significant speedups over non-preconditioned solvers for complex geometries and meshes specifically adapted to capture anisotropic features of a solution, including discontinuities arising from corners and edges.

Introduction

This paper is concerned with the efficient numerical solution of 3D wave propagation problems in largescale unbounded domains. The accurate numerical modeling of acoustic, elastodynamic or electromagnetic problems is a challenging task. It is also a timely research field due to the variety of possible applications, e.g., seismology, non-destructive testing, noise reduction, radar stealth, etc. We focus on an exterior Helmholtz problem for clarity, but the extension to other wave propagation problems (e.g., elastodynamics and electromagnetism) is straightforward. Let Ω ⊂ R 3 be a closed bounded domain with boundary Γ and let Ω + := R 3 \Ω denote the exterior scattering domain. Given an incident plane wave u i (x) = e ikd•x , x ∈ R 3 of wavelength λ (wavenumber k := 2π/λ) and direction d, the total field u = u i + u s (where u s is the scattered field) in the exterior domain satisfies the time-harmonic acoustic wave equation given by

∇ 2 u + k 2 u = 0 in Ω + . (1a)
We impose a sound-soft boundary condition as

u = 0 on Γ, (1b)
together with a Sommerfeld radiation condition on the scattered field u s = u -u i given as lim r→∞ ∂u s ∂r -iku s = 0, r = |x|.

Numerical methodologies to solve this system fall into two main categories: those obtained by domain methods requiring discretization of the entire physical domain of interest (e.g., finite difference [START_REF] Graves | Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences[END_REF], finite element [START_REF] Cohen | Finite element and discontinuous Galerkin methods for transient wave equations[END_REF][START_REF] Jin | The finite element method in electromagnetics[END_REF][START_REF] Khajah | Isogeometric finite element analysis of time-harmonic exterior acoustic scattering problems[END_REF][START_REF] Ganesh | High-order FEM domain decomposition models for high-frequency wave propagation in heterogeneous media[END_REF], finite volume and spectral methods [START_REF] Komatitsch | Introduction to the spectral element method for three-dimensional seismic wave propagation[END_REF]), and boundary methods [START_REF] Kress | Boundary integral equations in time-harmonic acoustic scattering[END_REF][START_REF] Liu | Recent advances and emerging applications of the boundary element method[END_REF][START_REF] Gaul | Boundary Element Methods for the Dynamic Analysis of Elastic, Viscoelastic, and Piezoelectric Solids[END_REF] requiring discretization of a reduced dimension such as the obstacle surface for an exterior scattering problem. Domain methods require a discretization of the exterior domain Ω + that leads to large but sparse algebraic systems. These systems are not well-suited to unbounded domains since the computational domain must be artificially truncated, leading to additional computational difficulties [START_REF] Clayton | Absorbing boundary conditions for acoustic and elastic wave equations[END_REF][START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Mattesi | A high-order absorbing boundary condition for 2D time-harmonic elastodynamic scattering problems[END_REF][START_REF] Modave | Corner treatment for high-order local absorbing boundary conditions in high-frequency acoustic scattering[END_REF]. In addition, the radiation condition (1c) must be approximated. A more suitable alternative is to recast the system (1) as a boundary integral equation (BIE). This subsequently leads to discretizations only on Γ which can be provided by a boundary element method (BEM). Where applicable, these methods are very efficient due to exact representations of both the boundary as well as radiation conditions at infinity. For a comparison between boundary element and finite element methods see, e.g., [START_REF] Harari | A cost comparison of boundary element and finite element methods for problems of time-harmonic acoustics[END_REF].

There exist various boundary integral equations to model acoustic wave scattering, with different advantages and drawbacks. We consider in this work the simplest approach. We will show that our methodology is nevertheless independent of the BIE considered. The solution u for an obstacle surface Γ can be sought as a boundary integral representation in the form of a single-layer potential operator S given by [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF] u(x) = Sq(x) = Γ Φ k (x, y)q(y)ds(y), x / ∈ Γ,

x ∈ Ω + , (2a)
where the fundamental solution of the 3D Helmholtz equation is given by

Φ k (x, y) = e ik|x-y| 4π|x -y| , x = y, k ∈ R + .
Expressing the unknown density q as q = -∂u/∂n (or, for geometries of no thickness, the jump [∂u/∂n] Γ across Γ) in the direction of the outer unit normal vector n to Γ, one can derive [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF][START_REF] Stephan | Boundary integral equations for screen problems in R 3[END_REF] a corresponding BIE for (2a) as

-u i = Γ Φ k (x, y)q(y)ds(y), x ∈ Γ. ( 2b
)
Clearly many realistic applications in acoustics (and other wave propagation problems) may also require Neumann or Robin boundary conditions. The methodology presented in this work can be straightforwardly extended to these cases. For the sake of simplicity, and without lose of generality, we restrict our simulations to the Dirichlet case.

A number of factors determine the effectiveness of a BEM applied for the solution of a general boundary integral equation. These include the choice of integral equation, their subsequent numerical resolution and the adopted discretization strategy [START_REF] Pérez-Arancibia | High-order integral equation methods for problems of scattering by bumps and cavities on half-planes[END_REF]. Standard implementations of BEM lead to dense and (possibly) non-symmetric linear systems whose solutions become prohibitively expensive for large-scale problems. This is particularly true for obstacles that are large relative to the wavelength and for problems employing highly non-uniform discretizations. The storage of such systems is of the order O(N 2 ), where N is the number of degrees of freedom corresponding to the discretization of the obstacle boundary (e.g., the number of nodes in the mesh). The computational complexity of solving such a dense system using a direct method like a Gaussian elimination is O(N 3 ), whereas resolution via an iterative method like GMRES is O(N iter N 2 ) if N iter is the number of iterations. Many approaches have been proposed to speedup the iterative resolution of these dense systems [START_REF] Ganesh | A spectrally accurate algorithm for electromagnetic scattering in three dimensions[END_REF][START_REF] Ganesh | A fast high order algorithm for multiple scattering from large sound-hard three dimensional configurations[END_REF][START_REF] Bruno | A fast, high-order algorithm for the solution of surface scattering problems: basic implementation, tests, and applications[END_REF][START_REF] Bruno | Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high frequency: the convex case[END_REF], among which Fast Multipole Methods (FMMs) have enjoyed considerable success in mechanical engineering problems [START_REF] Chaillat | A multi-level fast multipole BEM for 3-D elastodynamics in the frequency domain[END_REF][START_REF] Takahashi | A wideband fast multipole accelerated boundary integral equation method for timeharmonic elastodynamics in two dimensions[END_REF][START_REF] Cao | A fast directional BEM for large-scale acoustic problems based on the Burton-Miller formulation[END_REF][START_REF] Liu | Fast multipole boundary element method: theory and applications in engineering[END_REF] by enabling a fast evaluation of the matrix-vector product required by the iterative solver. Initially developed for N-body simulations, the FMM has since been extended to oscillatory kernels [START_REF] Darve | The fast multipole method: Numerical implementation[END_REF][START_REF] Greengard | Accelerating fast multipole methods for the helmholtz equation at low frequencies[END_REF][START_REF] Messner | Fast directional multilevel summation for oscillatory kernels based on Chebyshev interpolation[END_REF] and thus expanded its efficacy to many applications. An algebraic alternative designed for dense systems is based on the concept of hierarchical matrices (H-matrices) [START_REF] Bebendorf | Hierarchical Matrices: A Means to Efficiently Solve Elliptic Boundary Value Problems[END_REF]. The principle of H-matrices is to partition the initial dense linear system and reduce it to a data-sparse one by finding sub-blocks in the matrix that can be accurately approximated by low-rank matrices. The efficiency of hierarchical matrices relies on the possibility to approximate, under certain conditions, the underlying kernel function by low-rank matrices. The approach has been shown to be very efficient for asymptotically smooth kernels (e.g., the Laplace kernel) and relatively efficient in a large frequency range for oscillatory kernels such as Helmholtz or elastodynamic kernels [START_REF] Coulier | Application of hierarchical matrices to boundary element methods for elastodynamics based on Green's functions for a horizontally layered halfspace[END_REF][START_REF] Chaillat | Theory and implementation of H-matrix based iterative and direct solvers for Helmholtz and elastodynamic oscillatory kernels[END_REF][START_REF] Haider | Adaptive Cross Approximation for Bem in Elasticity[END_REF]. These two approaches have advantages and drawbacks in the context of the acceleration of BEMs. On the one hand, H-matrix based BEMs (H-BEM) are very easy to implement and have the advantage to be an algebraic approach. Most of the computational time of H-BEMs is spent to evaluate the data-sparse approximation of the BEM matrix, the cost of a matrix-vector is then very cheap. On the other hand, Fast Multipole accelerated BEMs (FM-BEMs) have been shown to have an optimal complexity in terms of computational times and memory requirements for wave propagation problems. The counterpart is a more involved implementation effort. In the following, our fast BEM solver is based on the FM-BEM for wave propagation problems.

Appropriate mesh adaptations can further improve the computational efficiency of the FM-BEM (particularly for complex geometries) by accommodating anisotropic features of a solution (e.g., some elastic materials) as well as discontinuities near geometric singularities (e.g., corners and edges). These behaviors are difficult to capture and ultimately diminish the regularity of the boundary solution and subsequent performance of a BEM. The underlying principle of a refinement algorithm relies on a transformation of an initial mesh into an improved one according to error estimates calculated at each step. The ultimate goal is a reduction in the number of degrees of freedom required to resolve a solution within a desired level of accuracy. Adaptation is particularly important for scattering obstacles that contain geometric singularities leading to rapid variations of surface solutions. If an extensive literature exists for volume-based methods [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF][START_REF] Gorynina | An easily computable error estimator in space and time for the wave equation[END_REF], relatively little research attention has been paid to error estimators and mesh adaptation techniques in the BEM community. In general, most adaptation techniques for BEMs have centered on isotropic techniques [START_REF] Feischl | Quasi-optimal convergence rate for an adaptive boundary element method[END_REF][START_REF] Bakry | A new accurate residual-based a posteriori error indicator for the BEM in 2D-acoustics[END_REF][START_REF] Bespalov | Adaptive BEM with optimal convergence rates for the Helmholtz equation[END_REF][START_REF] Kita | Error estimation and adaptive mesh refinement in boundary element method, an overview[END_REF] that improve performance but do not always recover optimal convergence rates for 3D problems with anisotropic features [START_REF] Aurada | Estimator reduction and convergence of adaptive BEM[END_REF]. In this work, we employ a recently introduced anisotropic mesh adaptation (AMA) strategy using a metric-based error estimator whose effectiveness in recovering convergence order was first demonstrated for volumetric (finite element) methods [START_REF] Loseille | Continuous mesh framework part I: well-posed continuous interpolation error[END_REF][START_REF]Continuous mesh framework part II: validations and applications[END_REF] and only recently for first-order BEMs [START_REF] Chaillat | Metric-based anisotropic mesh adaptation for 3D acoustic boundary element methods[END_REF]. The methodology is independent of discretization technique as well as the choice of partial differential equation (PDE) and integral equation formulation, iteratively constructing meshes refined in size, shape and orientation according to an "optimal" metric relying on a reconstructed Hessian of the boundary solution. The resulting adaptation is anisotropic in nature and accommodates geometric complexities that include engineering detail.

Since we are interested in large-scale scattering problems requiring iterative solutions, it is important to note that anisotropically (and even isotropically) adapted meshes worsen the conditioning of the overall FM-BEM system [START_REF] Graham | Anisotropic mesh refinement: the conditioning of Galerkin boundary element matrices and simple preconditioners[END_REF]. This leads to an increase in the number of iterations and thus overall computational time required to resolve the system. Even with the accelerations provided by the FMM, global solution times can still be prohibitively expensive, necessitating the use of preconditioning techniques to improve convergence of an applied iterative method. This is already an active area of research for uniform meshes [START_REF] Chaillat | A preconditioned 3-D multi-region fast multipole solver for seismic wave propagation in complex geometries[END_REF][START_REF] Takahashi | Application of the inverse fast multipole method as a preconditioner in a 3D Helmholtz boundary element method[END_REF]. We are interested here primarily by algebraic preconditioners which can be generally applied to a wide variety of problems as opposed to analytic preconditioners, which tend to be very efficient but problem specific [START_REF] Levadoux | A well-conditioned boundary integral equation for transmission problems of electromagnetism[END_REF][START_REF] Antoine | Analytic preconditioners for the boundary integral solution of the scattering of acoustic waves by open surfaces[END_REF][START_REF] Darbas | Combining analytic preconditioner and fast multipole method for the 3-D Helmholtz equation[END_REF][START_REF] Niino | Preconditioning based on Calderon's formulae for periodic fast multipole methods for Helmholtz' equation[END_REF]. This is motivated by the fact that FM-BEMs have been successfully applied to a wide variety of wave propagation problems, and the fact that the anisotropic mesh adaptation procedure used herein is independent of PDE and integral representation. Traditional algebraic preconditioning approaches such as incomplete LU, Sparse Approximative Inverse [START_REF] Carpentieri | Sparse pattern selection strategies for robust frobeniusnorm minimization preconditioners in electromagnetism[END_REF][START_REF] Carpentieri | Combining fast multipole techniques and an approximate inverse preconditioner for large electromagnetism calculations[END_REF], multi-grid methods [START_REF] Carpentieri | A matrix-free two-grid preconditioner for solving boundary integral equations in electromagnetism[END_REF] and nested GMRES algorithms [START_REF] Chaillat | A preconditioned 3-D multi-region fast multipole solver for seismic wave propagation in complex geometries[END_REF] have been applied to acoustic, electromagnetic or elastodynamic FM-BEMs. Since to reduce memory requirements, the FM-BEM does not provide access to the full system but only close interactions, these algebraic preconditioners have shown moderate efficiency. They probably do not contain enough information on the physics of the underlying continuous operator.

Our objective in this paper is to maintain the optimal complexity of the FM-BEM in terms of memory requirements and computational costs while complementing it with an algebraic preconditioner constructed from the complete BEM system. We are particulary interested in non-uniform meshes constructed by an anisotropic mesh adaptation procedure that is independent of the PDE and integral representation. In this context, our approach effectively utilizes the complementarities of the FM-BEM and the H-BEM to treat geometries modeled by non-uniform meshes. It is based on a nested sequence of GMRES iterations. Our methodology is presented as follows: in Section 2 we summarize the Boundary Element Method and its multi-level Fast Multipole-acceleration; in Section 3 we present the nested GMRES implementation using hierarchical matrices and consider some verification examples; and, finally, we demonstrate in Section 4 the efficacy of the iterative solver applied to anisotropic meshes through a variety of numerical experiments on an infinitely thin planar screen and a more complex F-15 aircraft.

Adaptive fast multipole-accelerated boundary element methods

We begin by presenting an adaptive FM-BEM for 3D wave propagation problems. Sections 2.1 and 2.2 summarize the general-order triangular Boundary Element formulation and subsequent multi-level Fast Multipole acceleration. Section 2.3 summarizes the mesh adaptation algorithm utilized to generate anisotropic meshes tailored to solutions of scattering problems on complex domains.

Boundary element methods

The numerical solution of the boundary integral equation (2b), known as the Boundary Element Method [START_REF] Bonnet | Boundary integral equations methods in solids and fluids[END_REF], is obtained by discretization of the surface Γ into N E finite (boundary) elements. Importantly in the BEM, each element E e must be analytically described. Each physical element E e on the approximate boundary is mapped onto a reference element ∆ e in the (ξ 1 , ξ 2 )-plane via

ξ ∈ ∆ e → y(ξ) ∈ E e , 1 ≤ e ≤ N e .
The boundary of the geometry is parametrically represented by using a number of d nodes and shape functions v e k , per element. We use as nodes y k all the vertices for a P 0 or P 1 mesh (d = 3) and all vertices plus midpoints for a P 2 mesh (d = 6) such that the geometry is approximated by

y(ξ) ≈ d k=1 y k v e k (ξ).
The shape functions are related to the canonical basis (v k ) 1≤k≤d defined on the reference element ∆ e by v e k (y(ξ)) = vk (ξ). Then, each unknown field q is also approximated on the element E e with interpolation functions (w i (y)) 1≤i≤N :

w i (y j ) = δ ij for 1 ≤ i, j ≤ N by q(y) ≈ d k=1 q k w e k (y),
where q k denotes the approximation of the nodal value of q(y k ). For isoparametric elements such as P 1 and P 2 , d = d , the shape functions serve as interpolation functions as well, and the N interpolation points coincide with the nodes previously defined. For a P 0 mesh, d = 1, the N interpolation points correspond to the centers of elements and we use piecewise constant interpolation functions.

To discretize the boundary integral equation (2b), we employ the collocation approach. It consists in enforcing the equation at a finite number of collocation points x. To have a solvable discrete problem, one has to choose N collocation points. The N interpolation points thus defined also serve as collocation points, i.e., (x i ) 1≤i≤N = (y j ) 1≤j≤N . In addition, a standard Gaussian quadrature formula is used to evaluate the integral in (2b). This discretization process transforms (2b) into a square complex-valued linear system of size N of the form Aq = b,

where the (N )-vector q collects the degrees of freedom (DOFs) while the (N )-vector b arises from the imposed incident wave field. Assembling the full dense matrix A requires the computation of all element integrals for each collocation point, thus requiring an O(N 2 ) computational time and memory. If an iterative method is used, each GMRES iteration requires one evaluation of Aq for a given q, a task requiring a computing time of order O(N 2 ). To lower this O(N 2 ) complexity, prohibitive for large BEM models, fast BEM solution techniques such as the Fast Multipole Method (FMM) must be employed.

Fast multipole-accelerated boundary element methods for wave propagation

The goal of the FMM is to accelerate the evaluation of the matrix-vector product Aq for a given q, required at each iteration of an iterative solver applied to the BEM-discretized equations. Moreover, the governing BEM matrix is never explicitly formed, which leads to a storage requirement well below O(N 2 ). Substantial savings in both computational time and memory requirements are thus achieved. In general terms, the FMM exploits a reformulation of the fundamental solution in terms of products of functions of x and of y, so that (unlike in the traditional BEM) integrations with respect to y can be reused when the collocation point x is changed. On decomposing the position vector r = yx into r = (y -y 0 ) + r 0 -(x -x 0 ), where x 0 and y 0 are two poles and r 0 = y 0 -x 0 , and invoking a well-known plane wave decomposition, the Helmholtz fundamental solution is written as [START_REF] Darve | The fast multipole method: Numerical implementation[END_REF]:

Φ k (x, y) = lim L→+∞ ŝ∈S e ikŝ.(y-y0) G L (ŝ; r 0 ; k)e -ikŝ.(x-x0) dŝ, ( 4
)
where S is the unit sphere of R 3 and the transfer function G L (ŝ; r 0 ; k) is defined in terms of Legendre polynomials P p and spherical Hankel functions of the first kind h

p by: (1) p (k|r 0 |)P p cos(ŝ, r 0 ) .

G L (ŝ; r 0 ; k) = ik 16π 2 0≤p≤L (2p + 1)i p h
It can be shown that expression ( 4) is valid only for well-separated sets of collocation and integration points clustered around poles x 0 and y 0 [START_REF] Darve | The fast multipole method: Numerical implementation[END_REF].

In the single-level FMM, a 3D cubic grid of linear spacing d embedding the whole boundary Γ is then introduced to drive the computation. The FMM basically consists in using decomposition (4), with the poles x 0 and y 0 being chosen as the cell centres of the cells C x and C y , whenever x and y belong to non-adjacent cubic cells. The treatment of such Fast Multipole (FM) contributions exploits the plane wave decomposition (4) of the fundamental solution, truncated at a finite L and in a manner suggested by its multiplicative form. When x and y belong to adjacent cells, traditional BEM evaluation methods are employed. To improve further the computational efficiency of the FM-BEM, standard (i.e., non-FMM) calculations must be confined to the smallest possible spatial regions while retaining the advantage of clustering the computation of influence terms into non-adjacent large groups whenever possible. This is achieved by recursively subdividing cubic cells into eight smaller cubic cells. New pairs of non-adjacent smaller cells, to which plane wave expansions are applicable, are thus obtained from the subdivision of pairs of adjacent cells. The cellsubdivision approach is systematized by means of an octree structure of cells. At each level , the linear cell size is denoted d . The level = 0, composed of only one cubic cell containing the whole surface Γ, is the tree root. The subdivision process is further repeated until the finest level = ¯ , implicitly defined by a preset subdivision-stopping criterion (d ¯ ≥ d min ), is reached. Level-¯ cells are usually termed leaf cells. This is the essence of the multi-level FMM, whose theoretical complexity is O(N log N ) per GMRES iteration both for computational time and memory requirements. More details on the implementation of the FM-BEM for waves can be found, e.g., in [START_REF] Chaillat | A multi-level fast multipole BEM for 3-D elastodynamics in the frequency domain[END_REF] where the method has been successfully employed in elastodynamics.

Anisotropic mesh adaptation for fast multipole-accelerated boundary element methods

The influence of the mesh discretization on the accuracy of a numerical solution still poses a challenge to the BEM community [START_REF] Bespalov | Adaptive BEM with optimal convergence rates for the Helmholtz equation[END_REF][START_REF] Bakry | A new accurate residual-based a posteriori error indicator for the BEM in 2D-acoustics[END_REF][START_REF] Chaillat | Metric-based anisotropic mesh adaptation for 3D acoustic boundary element methods[END_REF]. Anisotropic features of a solution (e.g., some elastic materials) as well as discontinuities near geometric singularities (e.g., corners and edges) are difficult to capture and ultimately diminish the regularity of the boundary solution and subsequent performance of a BEM. This is particularly true when uniform meshes are employed as singularities will degrade the convergence order. To this end, iterative mesh refinement schemes have been proposed to transform an initial mesh into an improved one according to error estimates calculated at each step, with the goal of reducing the number of degrees of freedom required to resolve a solution within a desired level of accuracy. Most work to date for BEMs has been confined to isotropic techniques based on a posteriori error analysis from which error indicators can be derived, see e.g., [START_REF] Bakry | A new accurate residual-based a posteriori error indicator for the BEM in 2D-acoustics[END_REF][START_REF] Bespalov | Adaptive BEM with optimal convergence rates for the Helmholtz equation[END_REF] for wave propagation problems. These indicators steer refinement by systematically marking and subdividing only elements where the error is above a specified threshold-a process known as Dörfler marking. Not only is this process unable to produce anisotropic discretizations, it is also heavily reliant on discretization strategy and model [START_REF] Erath | Simple error estimators for the galerkin BEM for some hypersingular integral equation in 2D[END_REF] (e.g., often requiring the more computationally expensive Galerkin approach). Finally, this strategy is not optimal for high-order discretizations [START_REF] Bespalov | Adaptive BEM with optimal convergence rates for the Helmholtz equation[END_REF].

We employ here a truly anisotropic metric-based boundary mesh adaptation recently introduced [START_REF] Chaillat | Metric-based anisotropic mesh adaptation for 3D acoustic boundary element methods[END_REF]. It does not use a Dörfler marking but rather generates a sequence of non-nested meshes. This is a nonlinear process that is solved iteratively. At each iteration, it seeks a new mesh that minimizes the L 2 linear interpolation error on the exact solution u.

A convenient framework to generate anisotropic meshes is that of Riemannian metric spaces. A Riemannian metric space is defined by a metric tensor M(x), i.e., is a smoothly varying function of the physical variable x. A surface mesh, associated with a Riemannian metric space M = (M(x)) x∈Γ , is a triangulation of the surface. Each mesh vertex x has an assigned value M(x) which dictates the size and orientation of adjacent elements. By generating a unit mesh in the corresponding Riemannian metric space, we obtain an anisotropic mesh refined in Euclidean space according to the metric M. This is the fundamental idea of metric-based mesh adaptation. Since the existence of a conforming unit mesh in which each triangle is perfectly unit with respect to a given Riemannian metric space is not guaranteed in general, the constraint must be relaxed somewhat. In practice, we seek a quasi-unit mesh with respect to M which is a mesh composed of quasi-unit triangles. A triangle, which is defined by its edges

{e i } i=1,3 , is quasi-unit with respect to M if 1 √ 2 ||e i || M , 1 ≤ i ≤ 3 and |K| M = √ 3 4 .
The principle of our adaptive mesh strategy is then to minimize the error coming from the linear interpolation of the (unknown) exact solution. The discrete optimization problem reads:

Find the surface mesh T opt N with N nodes, minimizing the discrete interpolation error ||u -

Π h u|| L p (Γ h ) , (6)
where Γ h is the discrete support representing Γ and Π h u is the linear interpolant of u on the mesh T opt N . This problem is generally intractable practically and would require some simplifications. The originality of our approach is to replace the discrete interpolation error with its continuous equivalent. Using the mathematical tools defined on the continuous mesh space instead of the discrete one, the optimization problem is simplified. In the continuous mesh setting [START_REF]Continuous mesh framework part II: validations and applications[END_REF], we need to introduce the complexity N of a mesh. It is the continuous counterpart of the number of vertices and it governs the size of the meshes, i.e.,

N := Γ d(x)dx, where d = (h 1 h 2 ) -1 is the density and (h i ) i=1,
We have shown [START_REF] Chaillat | Metric-based anisotropic mesh adaptation for 3D acoustic boundary element methods[END_REF] that the solution to [START_REF] Bebendorf | Hierarchical Matrices: A Means to Efficiently Solve Elliptic Boundary Value Problems[END_REF] is given by

M = N Γ det(|R s (H)|) 1 3 -1 det(|R s (H)|) -1 6 |R s (H)|, (8)
where R s is a restriction operator for surfaces and H is the Hessian of the solution. The mesh generated by this metric minimizes the L 2 -interpolation error.

The above results are valid for the linear interpolant Π h u of the exact solution u. However, neither Π h u or u are known. Only the approximation of u obtained by the FM-BEM is at hand. In order to relate the interpolation error to the approximation error, we consider a reconstruction operator R h which is applied to the numerical approximation u h . This operator can be a recovery process, a hierarchical basis, or an operator connected to an a posteriori estimate. There are numerous techniques in the literature to approximate or "recover" the Hessian (see, e.g., [START_REF] Picasso | A numerical study of some Hessian recovery techniques on isotropic and anisotropic meshes[END_REF][START_REF] Vallet | Numerical comparison of some Hessian recovery techniques[END_REF]). From a practical point of view, we only need to recover a pointwise 3D Hessian at the nodes of the current surface mesh. The method we employ here is based on a local quadratic recovery and uses the boundary solution only, see [START_REF] Chaillat | Metric-based anisotropic mesh adaptation for 3D acoustic boundary element methods[END_REF] for more details.

With the minimizing mesh in hand, one repeats the process for larger complexities until some ad-hoc stopping criterion is achieved. More precisely, given a user-prescribed sequence of complexity N = [N 1 , ..., N k ], we seek the sequence of corresponding optimal meshes T = [T 1 , ..., T k ]. This process is non-linear, i.e., both the mesh and the solution have to be converged. The following iterative algorithm is used to generate T i :

Start from the mesh T 0 i = T 0 i-1 (or from the initial uniform mesh T 0 at iteration 0). for all k = 0 : n iter do Compute the approximation u h on the mesh T k i with the iterative FM-BEM solver (software COFFEE). Compute the recovered Hessian from u h and deduce Metric [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF], with complexity N i . From [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF], a new quasi-unit mesh T k+1 i of complexity N i is generated with Feflo.a [START_REF] Loseille | Handbook of Numerical Methods for Hyperbolic Problems -Applied and Modern Issues[END_REF].

end for

The proposed strategy constructs adapted meshes that have been shown to recover optimal convergence rates for domains with corners and edges [START_REF] Chaillat | Metric-based anisotropic mesh adaptation for 3D acoustic boundary element methods[END_REF]. The main advantage of this strategy is to provide manageable system sizes for a desired accuracy (to achieve the same with a uniform mesh, one may require a discretization exceeding memory capabilities) to the FM-BEM solver. Importantly, the computation of the metric and subsequent generation of the mesh itself is a rather inexpensive calculation, negligible relative to the corresponding FM-BEM solution time (usually less than 5% of the overall computational cost). It is additionally independent of the PDE and integral equation formulation and hence can be applied to a wide variety of wave propagation problems. Figure 1 illustrates a corner portion of the three iterations of meshes generated for the planar screen geometry considered in Section 4.1. One can observe clear refinement at the edges. Unfortunately, such heavily non-uniform discretizations (note the stretched elements next to regular elements) lead to poorer conditioning of the FM-BEM system [START_REF] Graham | Anisotropic mesh refinement: the conditioning of Galerkin boundary element matrices and simple preconditioners[END_REF], further motivating the need for an efficient iterative solver.

Principles of the H-matrix representation of the BEM system for wave propagation

Hierarchical matrices or H-matrices have been introduced by Hackbusch [START_REF] Hackbusch | A Sparse Matrix Arithmetic Based on H-Matrices. Part I: Introduction to H-Matrices[END_REF] to compute a data-sparse representation of some special dense matrices (e.g., matrices resulting from the discretization of non-local operators). The principle of H-matrices is (1) to partition the matrix into blocks and (2) to perform low-rank approximations of the blocks of the matrix which are known a priori (via an admissibility condition) to be accurately approximated by low-rank decompositions. With these two ingredients it is possible to define fast iterative and direct solvers for matrices having a hierarchical representation. Using low-rank representations, the memory requirements and costs of a matrix-vector product are reduced.

The key component of hierarchical matrices is the recursive block subdivision. The first step prior to the partition of the matrix is a partitioning based on the geometry of the set of row and column indices of the system A. The objective is to permute the indices in the matrix in order to reflect the physical distance and hence interaction between degrees of freedom. Consecutive indicies should correspond to DOFs that interact at close range. A binary tree T I is used to drive the clustering where each node of the tree defines a subset of indices σ ⊂ I and each subset corresponds to a part in the partition of the domain. For each node we determine the box enclosing all the points in the cloud and subdivide it into 2 boxes, along the largest dimension. The subdivision is stopped when a minimum number of DOFs per box is reached (N leaf ).

After the clustering of the unknowns is performed, a block cluster representation T I×I of the matrix A is defined by going through the cluster tree T I . Each node of T I×I contains a pair (σ, τ ) of indices of T I and defines a block of A. This uniform partition defines a block structure of the matrix with a full pattern but is not optimal. As a matter of fact, some parts of the matrix A can be accurately approximated by a low-rank matrix at a higher level (i.e. for larger clusters). Such blocks are said to be admissible. A hierarchical representation P ⊂ T I×I that uses the cluster tree T I and the existence of admissible blocks is more appropriate. Starting from the initial matrix, each block is recursively subdivided until it is either admissible or the leaf level is achieved. For complex 3D geometries, an admissibility condition based on the geometry and the interaction distance between points is used to determine a priori the admissible blocks. In summary, the blocks of the partition can be of three types: at the leaf level a block can be either an admissible block or a non-admissible block, at a non-leaf level a block can be either an admissible block or an H-matrix (i.e a block that will be subsequently hierarchically subdivided).

H-matrix representations have been derived for some specific problems and will not result in efficient algorithms for all equations or matrices. The crucial point is to know a priori (1) if savings will be obtained when trying to approximate admissible blocks with a sum of separated variable functions and (2) which blocks are admissible since the explicit computation of the rank of all the blocks would be too expensive. In the case of asymptotically smooth kernels G(x, y), it is proved that under some a priori condition on the distance between x and y, the kernel is a degenerate function. More precisely, the admissibility condition is given by max

(diam X, diam Y ) ≤ η dist(X, Y ), (9)
where diam X (resp. diam Y ) is the diameter of a block domain

X (resp. Y ), dist(X, Y ) = inf{ x -y , x ∈ X, y ∈ Y }
is the Euclidean distance between the domains X and Y and η > 0 is the admissibility constant of the method (fixed to η = 3 in the following sections). After discretization, this property is reformulated as the efficient approximation of blocks of the matrix by low-rank matrices. The Laplace Green's function is an example of an asymptotically smooth kernel for which H-matrix representations have been shown to be very efficient. The Helmholtz Green's functions are not asymptotically smooth but oscillatory. It has been shown in [START_REF] Börm | Hierarchical matrices[END_REF] that even though the method is less optimal in that case, substantial memory savings are still achieved with a complexity level lower than O N 3/2 log N in practice.

Once the H-matrix representation of a matrix is computed, it is easy to derive an H-matrix based iterative solver. The only operation required is an efficient matrix-vector product. It is performed hierarchically by going through the block cluster tree. At the leaf level, there are two possibilities. If the block does not admit a low-rank approximation (non-admissible block), then the standard matrix-vector product is used. Otherwise, the block is marked as admissible such that a low-rank approximation (via adaptive cross approximation with stopping tolerance ε ACA corresponding to the accuracy of the low-rank representation) has been computed. The cost of this part of the matrix-vector product is then reduced. The interested reader is referred to [START_REF] Chaillat | Theory and implementation of H-matrix based iterative and direct solvers for Helmholtz and elastodynamic oscillatory kernels[END_REF] for full details concering the H-matrix construction for oscillatory Helmholtz kernels.

A preconditioned adapted FM-BEM solver

The objective of this work is to use the complementary aspects of the fast multipole method and H-matrix representations of the system. Our overall aim is to reduce total computational cost of the adaptive FM-BEM. Since the FM-BEM solvers have been shown to have an optimal complexity for wave propagation problems, our goal is to supplement the FM-BEM with a preconditioning provided by the H-matrix representation of the BEM system. The goal is to solve the right preconditioned FM-BEM system (3), i.e., find M such that

AM -1 Mq = b.
By doing so, we will solve AM -1 y = b and Mq = y with the goal of improving the conditioning of A. The best preconditioner is of course M = A but the cost of evaluating M -1 would be equivalent to solving the initial system. Instead, we want to use an approximation of A. In this work, we study the efficiency of the use of an H-matrix approximation of the BEM system to precondition the FM-BEM system in the context of anisotropic meshes. Since we only need a coarse approximation of M -1 , we will use a two level GMRES. The inner GMRES is used to calculate the preconditioner while the outer GMRES solves the corresponding preconditioned system. A flexible GMRES [START_REF] Saad | A flexible inner-outer preconditioned gmres algorithm[END_REF] allows the preconditioner to change at every iteration, hence permitting the use of nested "inner-outer iterations" which can invoke an iterative method itself as the preconditioner (e.g., another GMRES). Since the subsequent matrix-vector product is very cheap, an H-matrix based calculation is an ideal candidate for these inner iterations: it is algebraically formed and thus general, all the while respecting the integrity of the BEM system. Algorithm 1 summarizes the nested GMRES implementation employing such a preconditioner applied to system (3). The principle difference from the standard GMRES is found in line 3. Instead of the application of a fixed preconditioner M (e.g., right-preconditioning z j = M -1 v j ), we solve Aq j = v j approximately by a given iterative method-in this case, GMRES using (coarse) Hmatrix representations. This inner solve step preconditions the outer solve whose matrix-vector products are facilitated by the FMM.

Algorithm 1 Nested GMRES for preconditioning fast BEMs

1: r 0 = b -(Aq 0 ) FMM , β = r 0 2 =⇒ v 1 = r 0 /β
// evaluate the initial residual with FM-BEM 2: for j = 1, 2, . . . until convergence do // outer GMRES iterations with FM-BEM system 3:

Solve A Hmat z j := v j // inner GMRES with H-mat to define preconditioner 4:

w := (Az j ) FMM // FMM approximation 5:

for i = 1 to j do // Construction of the Krylov subspace 6:

h i,j := z i • w 7: w := w -h i,j z i 8:
end for 9:

h j+1,j = w 2 and v j+1 = w/h j+1,j 10:

y j = arg min y βe 1 -H(1 : j + 1, 1 : j)y 2 // Hij = hij, solve minimization pb. to find the update 11:

q j = q 0 + [z 1 , z 2 , . . . , z j ]y j // updated solution 12: end for The brief diagram of Figure 3 summarizes the full solver. Initial coarse meshes are constructed to have a density of approximately d λ = 5 collocation points per wavelength, leading to a small initial calculation of the boundary solution. The adaptive sequence computes the Hessian and corresponding metric, ultimately generating a quasi-unit mesh within this space. An FM-BEM is then applied to this new mesh to rapidly produce a solution for the next iteration. The metric calculation and mesh generation are very cheap-the complexity of our solver is dominated by the FM-BEM. Hence any improvement in computational cost is felt by every iteration of the adaptive mesh procedure.

Implementation details and validations of the preconditioner strategy

All the aforementioned algorithms have been implemented in the FORTRAN language and within the in-house accelerated boundary element software for elastodynamics and acoustics known as COFFEE. Numerical experiments were conducted on a high-performance computer equipped with four multicore Intel Xeon E7-4820 processors for a total of 32 cores. All anisotropic meshes in Section 4 are constructed by the metric-based adaptation procedure (Section 2.3) using Metrix and Feflo.a software packages from the Gamma3 team at INRIA Saclay [START_REF] Loseille | Handbook of Numerical Methods for Hyperbolic Problems -Applied and Modern Issues[END_REF][START_REF] Loseille | Continuous mesh framework part I: well-posed continuous interpolation error[END_REF][START_REF]Continuous mesh framework part II: validations and applications[END_REF].

As a verification of the FM-BEM code (whose FMM-speedups for piecewise linear elements have already been demonstrated in [START_REF] Chaillat | Metric-based anisotropic mesh adaptation for 3D acoustic boundary element methods[END_REF][START_REF] Chaillat | A multi-level fast multipole BEM for 3-D elastodynamics in the frequency domain[END_REF]), we consider diffraction by a sphere of radius R of an incident field with Input of an initial coarse mesh T 0 and solution Computation of optimal metric M i and generation of a new mesh T i

Step i

Resolution by preconditioned FM-BEM

for Each solver iteration do Computation of preconditioner with H-matrices Matrix-vector product via FMM on preconditioned vector end for

Resolution by FM-BEM

for Each solver iteration do Matrix-vector product via FMM end for direction d = (-1, 0, 0) T and a non-dimensional wavenumber kR = 10 for a sequence of uniformly-refined discretizations employing P 0 , P 1 and curved P 2 elements. The FM-BEM system is solved by using the simplest preconditioner, i.e., the right-diagonally preconditioned GMRES [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] algorithm which scales all the diagonal entries of the system matrix to 1. We hence solve the system (3) as

Output of the final adapted mesh and solution

AM -1 y = b, q = M -1 y ( 10
)
for a preconditioner M defined as

M ij = A ij , if i = j, 0, if i = j. ( 11
)
Figure 4 reports the L 2 -norm of the error, between the approximated surface density q N and the analytical density q for the integral equation on the boundary of a sphere for a fixed number of elements N E using P 0 , P 1 and P 2 interpolations and given by

∂Ω (q -q N ) 2 dx 1/2 ∂Ω q 2 dx 1/2 .
There is little error analysis for 3D collocation BEM solvers but there exist some empirical studies, see, e.g. [START_REF] Atkinson | The numerical solution of integral equations of the second kind[END_REF].

As this is a case of a smooth geometry, one readily appreciates the recovery of the expected behaviour of the L 2 convergence order in all three cases. We fix the GMRES tolerance to 10 -4 , which is a common choice in the BEM community. It has also been shown to produce the expected convergence rates in terms of errors between the discrete and exact physical solution, in problems of interest in this work [START_REF] Chaillat | Metric-based anisotropic mesh adaptation for 3D acoustic boundary element methods[END_REF]. Since [START_REF] Börm | Hierarchical matrices[END_REF] is a simple preconditioning technique commonly used for single-layer operators, we henceforth refer to solutions using this diagonalpreconditioner as the original ("regular") "FM-BEM" and to solutions using the preconditioning nested GMRES-implemented with H-matrix representations-as "preconditioned FM-BEM" (denoted by "prec FM-BEM" in plots). To verify the preconditioner implementation, we again consider the unit sphere with an incident field of direction d = (-1, 0, 0) T , but for non-dimensional wavenumbers kR = 17, 24, 30 and 43. We uniformly discretize the surface for each wavenumber using P 1 elements (with a density d λ of approximately 15 collocation points per wavelength) to obtain a corresponding series of meshes ranging in size from N ≈ 30 000 to N ≈ 180 000 vertices. Figure 5 presents an illustrative discretization as well as the solution (real part and complex amplitude) for the case of kR = 43. For each discretization, H-matrix representations are calculated with a coarse resolution of ε ACA = 10 -1 , admissibility condition ( 9) with η = 3 and a minimum block size stopping criterion N leaf = 100. These parameters are chosen empirically through a series of numerical experiments, a version of which is presented in this section for uniformly discretized planar screens and later for anisotropically discretized planar screens. Each case is iteratively resolved either by a (right) diagonal-preconditioned FM-BEM (11) ("FM-BEM") or the nested H-matrix preconditioned FM-BEM ("prec FM-BEM"), with the (outer) GMRES tolerance fixed to 10 -4 (500 maximum iterations) and the inner tolerance of the nested GMRES fixed to 10 -1 (10 maximum inner iterations). Defining the computational setup time for the FM-BEM system and H-matrix representations (for an approximation tolerance ε = ε ACA ) as t FMM setup and t Hmat,ε setup , respectively, the corresponding setup (precompu-tation) times are given for FM-BEM and preconditioned FM-BEM respectively as

t setup = t FMM setup , t prec setup = t FMM setup + t Hmat,ε setup . (12)
The left plot of Figure 6 presents the corresponding times for each discretization and non-dimensional wavenumber of the sphere example. The addition of an H-matrix construction for the nested preconditioner results in a slightly longer precomputation time. Although the complexity of constructing the H-matrix representation is generally higher than the FMM complexity, it can be shown [START_REF] Bebendorf | Hierarchical Matrices: A Means to Efficiently Solve Elliptic Boundary Value Problems[END_REF][START_REF] Chaillat | Theory and implementation of H-matrix based iterative and direct solvers for Helmholtz and elastodynamic oscillatory kernels[END_REF] that the growth in complexity is bounded as O(N 3/2 log N ). This property is also verified in the figure where we have additionally plotted the setup time t Hmat,ε setup . The right plot gives corresponding costs of a single (outer) iteration of FM-BEM and preconditioned FM-BEM defined as

t iter = t FMM iter , t prec iter = t FMM iter + N Hmat iter t Hmat iter , (13)
where N Hmat iter is the number of iterations of the inner sequence of GMRES facilitated by H-matrix representations. One similarly observes that the single iteration times are longer for the nested GMRES since they include FMM matrix-vector computations as well as inner H-matrix computations (which are themselves rather cheap and are applied a maximum of ten times an outer iteration). The maximum numerical rank among all admissible blocks in the H-matrix representation is provided in parentheses and grows with complexity as expected. We then report in Table 1 the global time-to-solution given for each case by

t global = N iter t FMM iter + t FMM setup , t prec global = N prec iter t FMM iter + N Hmat iter t Hmat iter + t FMM setup + t Hmat,ε setup . (14)
In spite of these marginal increases in setup and single iteration times, it is significantly less for the preconditioned FM-BEM owing to the significant reduction in the number of overall (outer GMRES) iterations N prec iter when resolving the FM-BEM system. Defining the corresponding speedup S of preconditioned FM-BEM versus FM-BEM,

S = t global -t prec global t global × 100,
one already observes during the verification the larger than 50% speedup in computation time when employing the preconditioned solver for uniformly discretized spheres. The computational gains are demonstrated in Section 4 to be even more significant for the anisotropically-discretized meshes and higher-order numerical quadratures that are of interest. Geometric singularities and the subsequent discontinuous behavior of the corresponding boundary solution can lead to a significant performance loss in a BEM. An infinitely thin planar screen presents a challenging example in this context: it has sharp and thin geometric features that encourage, like many realistic configurations, a mesh adaptivity suitable to capture the behavior of the wave propagation around these geometric singularities. This is a case of particular interest in Section 4.1 and so in this spirit we consider a series of verifying numerical examples for scattering of plane waves of direction d = (1, 0, 0) T by a uniformly discretized screen of characteristic length, i.e., edge length L. As in for the sphere, we consider the screen uniformly discretized by piecewise-linear elements as seen in Figure 7, which presents an illustrative discretization as well as a solution (real part and complex amplitude) for the case of kL = 43. A density d λ of approximately 15 collocation points per wavelength is chosen to obtain a corresponding series of meshes ranging in size from N ≈ 300 000 to N ≈ 500 000 vertices. The corresponding setup and single-iteration times are presented in Figure 8. Again we can observe that the addition of an H-matrix construction for the nested preconditioning results in a slightly longer precomputation time due to the H-matrix representation, but is coherent with theoretical complexity. Similarly again, the iteration times are longer as well since they include an FMM matrix-vector computation as well as inner H-matrix computations (the numerical rank of the largest approximated admissible block is included in parentheses). The global time-to-solution given in Table 2, however, is significantly less for the preconditioned FM-BEM solver owing to the significant reduction in the number of overall (outer) GMRES iterations N prec iter when resolving the FM-BEM system, enabling more than 50% reduction in the global computation time.

N

We now consider the efficacy of the nested GMRES-based preconditioner with respect to the parameters most important to its additional computational time, i.e., that of the construction of the H-matrix representation. In order to optimize performance, a balance must be struck between precomputational cost and the number of iterations during which the corresponding matrix-vector product is applied. Similar studies are presented for anisotropic discretizations of the screen in Section 4.1. For each uniform discretization we solve the resulting system via FM-BEM and the preconditioned FM-BEM system for various non-dimensional wavenumbers and various accuracies for the low-rank approximations in the H-matrix representation. For a fixed non-dimensional wavenumber kL = 50 and low-rank approximation tolerances of

ε ACA = 10 -1 , 2 • 10 -1 , 5 • 10 -1 ,
Figure 9 presents setup and single iteration times for a series of meshes ranging in size from N ≈ 30 000 to N ≈ 200 000 (these are no longer of fixed density). These parameters are of interest because they are very coarse approximations and are relatively quick to precompute. One can observe similar performance for each tolerance ε ACA , with the finest tolerance ε ACA having the shortest single iteration time of the resulting nested GMRES iterative solver. Fixing ε ACA = 10 -1 , we consider construction times and subsequent speedups for varying non-dimensional wavenumbers kL = [START_REF] Coulier | Application of hierarchical matrices to boundary element methods for elastodynamics based on Green's functions for a horizontally layered halfspace[END_REF][START_REF] Loseille | Handbook of Numerical Methods for Hyperbolic Problems -Applied and Modern Issues[END_REF]100].

It is known that the H-matrix representations in their current form are not optimized for very high nondimensional wavenumbers [START_REF] Chaillat | Theory and implementation of H-matrix based iterative and direct solvers for Helmholtz and elastodynamic oscillatory kernels[END_REF]. Table 3, which presents speedups over unpreconditioned FM-BEM, indicates however that this possibly is not the case when employing a coarse approximation only as a preconditioner; indeed, one can appreciate the significant speedups in all three cases, including for increasing wave numbers. Remark 3.1. The experiments above have employed modest numbers of Gaussian quadrature points to perform the numerical integrations of triangles required to construct the BEM system. For adapted meshes containing highly stretched, flat and anisotropic elements, this can be cumbersome as there may not be an adequate distribution of points, ultimately leading to poorer convergence of the iterative system. A larger number (13 points) of quadrature points may need to be used and are, in fact, required to maintain the properties of using higher-order element discretizations in a BEM. The rest of the experiments in this paper, unless otherwise noted, utilize a higher order of Gaussian numerical quadrature. This may increase certain computational costs but will reduce the number of iterations. Furthermore, as demonstrated in Table 4 for the above sphere and screen examples, the speedup obtained by the preconditioned FM-BEM may be even greater for more accurate quadrature routines.

N

Efficiency of the preconditioner in the presence of an anisotropic mesh

The sound-soft boundary condition (1b), when compared to other kinds of boundary conditions, gives rise to more severe singularities in the boundary solution when considering scattering by obstacles with corners or edges. This suggests the use of an anisotropic mesh adaptation precedure to account for the loss of regularity of the corresponding boundary solutions and hence order of convergence. Unfortunately, both isotropically and anisotropically adapted meshes worsen the conditioning of the overall BEM system, leading to an increase in the number of iterations and thus overall computational time required to iteratively resolve the solution [START_REF] Graham | Anisotropic mesh refinement: the conditioning of Galerkin boundary element matrices and simple preconditioners[END_REF]. This naturally motivates the use of an effective preconditioner to accommodate the subsequent increase in computational cost. Furthermore, our anisotropic mesh adaptation procedure is itself iterative and requires several FM-BEM solves; hence any decrease in computational time is felt multifold. In what follows, we consider the efficacy of the H-matrix based nested GMRES preconditioning on such adapted meshes constructed for problems of scattering by a planar screen and a more complex F-15 aircraft.

Diffraction by a screen

For the thin planar screen example of Section 3.3, the usage of an adaptation algorithm and, in particular an anisotropic one, is prudent: it is well known (see, e.g., [START_REF] Juhl | A note on the convergence of the direct collocation boundary element method[END_REF]) that at the edge of a screen, the jump behaves as ∂u ∂n ∼ (kr) -0.5 as kr → 0,

where r is the distance from the edge. At the corner of the screen, this singularity is more severe, and has been shown [START_REF] Bažant | Three-dimensional harmonic functions near termination or intersection of gradient singularity lines: a general numerical method[END_REF] to take the form ∂u ∂n ∼ (kr) -0.704 as kr → 0.

Owing to the singular behavior, the approximated solution on a sequence of uniformly refined meshes should yield a sub-optimal convergence rate. Hence planar screens are an excellent example of the capabilities of an (iterative) adapted mesh strategy, an example of which is given in Figure 1 for a corner of the screens used in the subsequent experiments. One observes in these images highly stretched and anisotropic elements near the corners and edges. These ultimately degrade the performance of an iterative solver, and hence it is of great interest to study preconditioning in this context.

ε ACA = 10 -1 , 2 • 10 -1 , 5 • 10 -1 ,
for a sequence of anisotropic meshes constructed with the mesh adaptation procedure for an incident plane wave of non-dimensional wavenumber kL = 10. There are expected slight increases in both setup time and iteration time for each choice of ε ACA . We observe each of these cases is almost identical in cost (since they are very similar in the maximum low rank approximation), and so we adopt henceforth ε ACA = 10 -1 in order to remain consistent with parameters found for uniform discretizations.

In addition, we report in Table 5 the characteristics of the anisotropic meshes, i.e., maximum and minimum length size together with anisotropic ratios. If we define γ K to be the ratio of the maximum edge length to minimum edge length in a triangle K, then we define the anisotropic ratio γ of a mesh T N as

γ max = max K∈T N γ K , γ mean = 1 #T N K∈T N γ K ,
where one can observe that the more anisotropic the mesh, the larger the value of γ K . It should be noted that these experiments are conducted on several meshes of a single adaptation sequence (i.e., as discretization increases so does the level of anisotropy). Figure 11 indicates setup times and single iteration times for the cases of kL = 50, 100 and 200, respectively. Each case carries a different series of meshes which are refined according to the mesh adaptation procedure for their respective wavelengths. Again, we observe slight increases in setup time, but the iteration time increases with level of anisotropy and size of the grid. Table 6 shows that, although with the more marginal increases in computation time, the method achieves even more significant speedups than that of a uniform discretization: on the order of 70% in most cases. This demonstrates the effectiveness of the preconditioned FM-BEM solver. One additionally observes a smaller speedup for finer (and more anisotropic discretizations). This is due to the non-optimality of the H-matrix representation for non-uniform discretizations. For example, the stopping criterion N leaf (minimal box size in the subdivision) must be tweaked to accommodate highly stretched elements. Finally, it is important to recall that we have used in all examples a simple diagonal-preconditioner even in the solutions denoted "FM-BEM". The solutions corresponding to the nested GMRES implemented with H-matrix representations are denoted "prec FM-BEM". In Table 7 we report also the number of iterations and solution times for the diffraction by an anisotropically discretized screen with kL = 50. These results illustrate the already important improvement of the conditioning due to the diagonal preconditioner, leading to an even better behaviour of our nested GMRES preconditioner.

ε = 10 -1 ε = 2 • 10 -1 ε = 5 • 10 -1 N Min length

N

Min

Conclusions and outlook

We have presented a nested GMRES-based adaptive fast BEM to model 3D wave propagation problems. It relies on an inner GMRES sequence of coarse and very fast H-matrix approximated matrix-vector products to precondition the Fast Multipole accelerated Boundary Element system used in the outer GMRES. We have demonstrated a reduction in the (outer) GMRES iterations and hence the subsequent invocations of the more expensive FM-based matrix-vector product, resulting in significant global time-to-solution speedups versus conventionally-preconditioned GMRES methods. In particular, we have studied the preconditioned solver applied to the problem of diffraction by an infinitely thin planar screen and an F-15 aircraft, both for uniform and highly non-uniform discretizations. The algorithm has been shown to provide significant speedups for most configurations including up to 80% for some anisotropic meshes.

Efficient construction and application of the H-matrix representations in general is still very much a work in progress for wave propagation. The current implementation is promising and is already useful to improve the cost of a mesh adaptation procedure. Further work will involve adapting the H-matrix construction specifically for highly anisotropic meshes including those with curved elements. An obvious first step would be to alter the distance function for partitioning based on the anisotropy of the mesh-this information can be provided by a metric that is derived from the adaptation procedure. And, finally, the methodology herein can be extended for double-and combined-layer formulations of the boundary integral equations, as well as extended for more complex physics such as elastodynamics problems; both the H-matrix approximation and the anisotropic mesh procedure are independent of choice of PDE and integral equation. Diffraction by an F-15 aircraft: number of solver iterations and global solution times (including setup) for kL = 114 using 3 Gaussian quadrature points in regular integrations over triangles. The number of iterations reported here is slightly larger than those reported in [START_REF] Chaillat | Metric-based anisotropic mesh adaptation for 3D acoustic boundary element methods[END_REF]; this is due to the latter use of an initial guess which is not employed in this experiment.

N

(Grant number 2015 60 0006).

  2 are the local sizes along the principal directions of the metric. The global optimization problem becomes to find the optimal continuous mesh M * minimizing the continuous interpolation error e Rs(M) in the L p norm: Given a complexity N , find M * = argmin M Γ e Rs(M) (x)

Figure 1 :

 1 Figure 1: An example of the first, third and fifth iterations (in the mesh adaptation process) employed in the anisotropic screen meshes considered in Section 4.1. The complexities input for each case were N = 350, 1400 and 5600, respectively.

  3. A nested GMRES-based preconditioner for adapted FM-BEM Section 3.1 briefly discusses the theory, construction and deployment of an H-matrix representation of the BEM system for wave problems. It is subsequently incorporated into a nested GMRES iterative solver described in Section 3.2. Section 3.3 provides several validating examples of the implementation of both the FM-BEM solver and the preconditioner constructed by H-matrix representations.

Figure 2 :

 2 Figure 2: The nested GMRES employing fast-multipole approximations preconditioned by a GMRES sequence of coarse Hmatrix approximations.

Figure 3 :

 3 Figure 3: Summary of the adaptive FM-BEM solver with preconditioning.

Figure 4 :

 4 Figure 4: Diffraction by a uniformly discretized sphere: convergence results of the computed density on the boundary for kR = 10 for discretizations with P 0 , P 1 and curved P 2 elements.

Figure 5 :

 5 Figure 5: Diffraction by a uniformly discretized sphere: a portion of the mesh; the real part of the boundary solution; and the complex amplitude of the boundary solution for kR = 43 (N ≈ 178 000).

Figure 6 :

 6 Figure 6: Diffraction by a uniformly discretized sphere: setup times (tsetup, t prec setup ) and single iteration times (t iter , t prec iter ) applied to meshes of fixed density d λ = 15. The maximum numerical rank along all admissible blocks in the H-matrix representation is given in parentheses.

Figure 7 :

 7 Figure 7: Diffraction by a uniformly discretized screen: a portion of the mesh; the real part of the boundary solution; and the complex amplitude of the boundary solution for kL = 43 (N ≈ 30 000).

Figure 8 :

 8 Figure 8: Diffraction by a uniformly discretized screen: setup times (tsetup, t prec setup ) and single iteration times (t iter , t prec iter ) applied to meshes of fixed density d λ = 15. The maximum numerical rank along all admissible blocks in the H-matrix representation is given in parentheses.

N

Figure 9 :

 9 Figure 9: Diffraction by a uniformly discretized screen: setup times (tsetup, t prec setup) and single iteration times (t iter , t prec iter ) for non-dimensional wavenumber kL = 50 and different H-matrix approximation tolerances ε ACA . The maximum numerical rank along all admissible blocks in the H-matrix representation is given in parentheses.

Figure 10 :

 10 Figure 10: Diffraction by an anisotropically discretized screen: setup times (tsetup, t prec setup) and single iteration times (t iter , t prec iter ) for non-dimensional wavenumber kL = 10 and different H-matrix approximation tolerances ε ACA . The maximum numerical rank along all admissible blocks in the H-matrix representation is given in parentheses.

Figure 10

 10 Figure 10 presents the setup time and single iteration times for three levels of accuracy in H-matrix representations, i.e. ε ACA = 10 -1 , 2 • 10 -1 , 5 • 10 -1 ,

Figure 11 :

 11 Figure 11: Diffraction by an anisotropically discretized screen: setup times (tsetup, t prec setup ) and single iteration times (t iter , t prec iter ) for non-dimensional wavenumbers kL = 50, 100 and 200. The maximum numerical rank along all admissible blocks in the Hmatrix representation is given in parentheses.

Figure 13 :

 13 Figure 13: Diffraction by an anisotropically discretized F-15 aircraft: setup times (tsetup, t prec setup) and single iteration times (t iter , t prec iter ) for non-dimensional wavenumber kL = 114. The maximum numerical rank along all admissible blocks in the H-matrix representation is given in parentheses.

Table 1 :

 1 Diffraction by a uniformly discretized sphere: number of solver iterations and global solution times (including setup) for meshes of fixed density d λ = 15.

		kR N iter	t global	N prec iter	t prec global	S
	27 870	17	84	9m56s	23	4m14s 57%
	54 117	24	89	25m32s	22 10m05s 61%
	87 615	30	115	59m00s	32 26m36s 55%
	178 472	43	209 188m09s	58 84m09s 55%

Table 2 :

 2 Diffraction by a uniformly discretized screen: number of solver iterations and global solution times (including setup) for meshes of fixed density d λ = 15.

Table 3 :

 3 Diffraction by a uniformly discretized screen: number of solver iterations and global solution times (including setup) for non-dimensional wavenumbers kL = 25, 50 and 100.

Table 4 :

 4 Number of solver iterations and global solution times (including setup) for discretizations of fixed density d λ = 15 for the sphere (top) and screen (bottom) using 13 Gaussian quadrature points in integrations over triangles.

Table 5 :

 5 Diffraction by an anisotropically discretized screen: number of solver iterations and speedups for non-dimensional wavenumber kL = 10 and varying accuracies ε ACA .

Table 6 :

 6 Diffraction by an anisotropically discretized screen: number of solver iterations and global solution times (including setup) for non-dimensional wavenumbers kL = 50, 100 and 200.

  .96 • 10 -3 9.51 • 10 -1 2.03 • 10 2 3.21 • 10 2 210 59m06s 18 10m19s 83% 15 687 1.97 • 10 -3 8.23 • 10 -1 3.22 • 10 2 4.18 • 10 2 180 66m48s 21 18m23s 73% 25 428 9.75 • 10 -4 1.14 • 10 -0 7.44 • 10 2 1.17 • 10 3

	length Max length	γ mean	γ max	N iter	t global	N prec iter	t prec global	S
	10 562 2180 97m41s	25 38m22s 61%

Table 8 :

 8 Diffraction by an F-15 aircraft: number of solver iterations and global solution times (including setup) for nondimensional wavenumber kL = 114.

  .96 • 10 -3 9.51 • 10 -1 2.03 • 10 2 3.21 • 10 2 395 59m54s 26 10m32s 82% 15 687 1.97 • 10 -3 8.23 • 10 -1 3.22 • 10 2 4.18 • 10 2 527 93m27s 32 19m26s 79% 25 428 9.75 • 10 -4 1.14 • 10 -0 7.44 • 10 2 1.17 • 10 3

	Min length Max length	γ mean	γ max	N iter	t global	N prec iter	t prec global	S
	10 562 2726 169m49s	76 79m27s 53%

Table 9 :

 9

Acknowledgements

This work has been funded by the Agence Nationale de la Recherche (RAFFINE project, ANR-12-MONU-0021). The research of F. Amlani has been partially funded by DGA, Direction Générale de l'Armement

Diffraction by an F-15 aircraft

We revisit the geometry of the unarmed F-15 aircraft (of characteristic length L in the longest direction) previously considered by the mesh adaptation algorithm in [START_REF] Chaillat | Metric-based anisotropic mesh adaptation for 3D acoustic boundary element methods[END_REF] for an incident plane wave with direction d = (1/2, 1, 0) T and non-dimensional wavenumber kL = 114. A series of adapted meshes are considered from a sequence of 15 that are generated (specifically, the third, sixth and ninth). Each mesh is an increasingly anisotropic discretization by piecewise-linear elements as seen in Figure 12, which presents an illustrative discretization of the anisotropic adaptation on a portion of the wing as well as the boundary solution (real part and complex amplitude). We see the high directionality of the diffracted waves with different patterns for the front and back of the aircraft, and most of the resulting anisotropy in the meshes are reached in areas where the diffracted waves have low regularity. The corresponding setup and single-iteration times are presented in Figure 13. Again we can observe that the addition of an H-matrix construction for the nested preconditioning results in a slightly longer precomputation time due to the H-matrix representation. Similarly again, the single iteration times are longer as well since they include an FMM matrix-vector computation as well as inner H-matrix computations. The largest approximated admissible low-rank block is included in parentheses. The global time-to-solution given in Table 8, however, is significantly less for the preconditioned FM-BEM solver owing to the significant reduction in the number of overall (outer) GMRES iterations when resolving the FM-BEM system, enabling up to more than 80% reduction in global computation time, significantly better performance than the uniform discretizations treated in Section 3.3.

Remark 4.1. As discussed in Remark 3.1, a number of 13 points was used to perform numerical integration for the experiments in this section. For the sake of consistency, we can also consider the parameters of [START_REF] Chaillat | Metric-based anisotropic mesh adaptation for 3D acoustic boundary element methods[END_REF] where three points are used for non-singular integration over a triangle. The resulting linear system requires significantly more iterations to resolve as seen in Table 9 due to the stretching of elements. Indeed, even in this case, speedups up to 80% are achieved when employing the preconditioned FM-BEM.