Philippe Marchner

Hadrien Bériot

Sophie Le Bras

Xavier Antoine

A domain decomposition solver

Keywords: Absorbing boundary conditions, optimized Schwarz method, high-performance computing, high frequency problems, Pierce operator, domain decomposition AMS subject classifications. 35J05, 65N55, 68W10

1. Introduction. The reduction of noise emissions in the transport industry has become of increasing importance in the past decades. Noise simulation tools are crucial to support experimental techniques and lower the design cost of quiet technologies. The present work is motivated by aeroacoustics applications. Solving all the scales of the aeroacoustic problem requires tremendous numerical resolution and modeling effort. We focus instead on a hybrid approach, where the mean flow and acoustic perturbations are solved in separated steps. Such an approach is relevant in the context of aeroacoustics analogies (see e.g. [START_REF] Lighthill | On sound generated aerodynamically I. General theory[END_REF][START_REF] Mohring | A well posed acoustic analogy based on a moving acoustic medium[END_REF]), in which acoustic sources are extracted from a mean flow computation in a first step, and are propagated thanks to a wave-like partial differential equation in a second step. Various wave propagation models exist depending on the mean flow assumptions and presupposed flow-acoustic interactions. One of the most complete model for non-viscous flows without heat conduction and linear acoustic perturbations is the linearized Euler equations (LEEs). The LEEs support acoustic, hydrodynamic and entropy waves and their linear interactions. It is a vector partial differential equation describing five acoustic variables in three-dimensions, that remains costly and is subjected to numerical instabilities. We focus our study on the Pierce operator [START_REF] Pierce | Wave equation for sound in fluids with unsteady inhomogeneous flow[END_REF], which has the advantage to be described by a single scalar variable. It can provide an accurate alternative to the LEEs, especially in the high frequency limit [START_REF] Bras | Direct and adjoint problems for sound propagation in non-uniform flows with lined and vibrating surfaces[END_REF]. Derived for arbitrary steady base flow, the operator is more general than the linearized potential equation from Goldstein [START_REF] Goldstein | Aeroacoustics[END_REF] which is limited to potential homentropic mean flows.

The Pierce operator has the advantage to be self-adjoint. Its resolution can be used as a black-box tool for optimization, compute tailored Green's functions or solve inversion problems such as sound source localization. We refer to the study initiated in [START_REF] Bras | Direct and adjoint problems for sound propagation in non-uniform flows with lined and vibrating surfaces[END_REF] for more details.

The Pierce operator takes the form of a convected wave operator with variable coefficients

P = -ρ 0 D Dt 1 ρ 2 0 c 2 0 D Dt + ∇ • 1 ρ 0 ∇ , D Dt = ∂ t + v 0 • ∇,
and v 0 (x) is the vector velocity field, which is assumed to be subsonic so that the condition ∥v 0 (x)∥/c 0 (x) < 1 holds. These quantities are data of the problem, which are computed in advance by a stationary flow solver, typically a Reynolds-averaged Navier-Stokes (RANS) solver. For a potential mean flow a non-linear Poisson solver can be used [START_REF] Gelder | Solution of the compressible flow equations[END_REF].

For turbomachinery applications it is customary to solve such a PDE in the frequency domain, because fixed frequency tones are often the main noise contribution [START_REF] Astley | Computational aero-acoustics for fan duct propagation and radiation. Current status and application to turbofan liner optimisation[END_REF]. Moreover porous material or acoustic liners are difficult to model in the timedomain. For all these reasons an efficient and accurate frequency-domain resolution of such an operator is critical, in particular at high frequencies.

At the mathematical level, time-harmonic wave propagation problems are notoriously difficult to solve at high frequencies [START_REF] Ernst | Why it is difficult to solve Helmholtz problems with classical iterative methods[END_REF]. The usual variational formulation of such problems are complex-valued and strongly indefinite. The strong indefinitness of the problem leads to inefficient iterative solvers, and makes preconditioners challenging to design [START_REF] Erlangga | Advances in iterative methods and preconditioners for the Helmholtz equation[END_REF]. As a result, direct solvers are generally advocated in the industry. High frequency problems face two major numerical difficulties: i) the number of discretization points per wavelength must grow significantly to avoid the pollution effect, especially for problem with flow convection [START_REF] Bériot | Analysis of high-order finite elements for convected wave propagation[END_REF]; ii) the memory requirement for direct solvers grows roughly as O(ω 3 ) for a fixed numerical resolution, making high frequency problems out of reach in practice. The first difficulty has been significantly reduced thanks to the progress of high-order discretization methods. In this work we use a conformal high-order finite element discretization based on integrated Legendre polynomials [START_REF] Solin | Higher-Order Finite Element Methods[END_REF]. Order adaptivity rules have been designed [START_REF] Bériot | Efficient implementation of high-order finite elements for Helmholtz problems[END_REF][START_REF] Bériot | Anisotropic adaptivity of the p-FEM for time-harmonic acoustic wave propagation[END_REF] to leverage the pollution effect. The counterpart of using a high-order basis is the increase of the matrix filling, which induces a higher memory footprint.

The second difficulty is common to all discretization methods, where the computer memory dictates the high frequency limit that can be simulated. Solutions to time-harmonic PDEs are non-local in space and highly oscillatory, making most of numerical algorithms hard to parallelize. Direct solvers can take advantage of the linear system structure and use efficient distributed parallelism strategies based on graph partitioning methods [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF]. While it provides significant speed ups, memory overhead and scalability limits the applicability of direct solvers to a large number of processes [START_REF] Li | 3D frequency-domain elastic wave modeling with the spectral element method using a massively parallel direct solver[END_REF]. This is expected since a pure algebraic parallel approach requires the factorization of a fully populated Schur complement system. Improvements of direct solvers were achieved in the recent years thanks to Block-Low Rank compression [START_REF] Amestoy | Performance and scalability of the block low-rank multifrontal factorization on multicore architectures[END_REF] combined with mixed-precision arithmetic [START_REF] Amestoy | Mixed precision low-rank approximations and their application to block lowrank LU factorization[END_REF]. The applicability to large scale time-harmonic problems is an active research topic.

The purpose of this article is to address the second difficulty based on a domain decomposition strategy at the continuous level. We reformulate the Pierce radiation problem using an optimized Schwarz domain decomposition method [START_REF] Gander | Optimized Schwarz methods without overlap for the Helmholtz equation[END_REF][START_REF] Antoine | Optimized Schwarz domain decomposition methods for scalar and vector Helmholtz equations[END_REF] posed on the partition interfaces, and propose a hybrid direct-iterative solver suitable for high frequency computations on a distributed memory architecture. The computational domain is split into independent non-overlapping subdomains where a direct solver is called, and the global interface problem is solved in parallel by an iterative method.

In that way we take advantage of direct solvers for medium size problems, and rely on solving iteratively an interface problem of smaller size. The resolution of the interface problem amounts to solve a Schur complement system at the algebraic level.

Domain decomposition also inherits from the strong indefiniteness of time-harmonic problems. The original Schwarz method does not converge for such problems and the issue was resolved in the early 1990's by Lions [START_REF] Lions | On the Schwarz alternating method. III: a variant for nonoverlapping subdomains[END_REF] and Després [START_REF] Després | Méthodes de Décomposition de Domaine pour les Problèmes de Propagation d'Ondes en Régime Harmonique[END_REF]. An appro-This manuscript is for review purposes only.

priate choice of transmission conditions between the subdomains is required for fast convergence, which depends directly on the underlying PDE. The ideal transmission condition is given by the Dirichlet-to-Neumann (DtN) map relative to the complementary of each subdomain, which is the continuous equivalent of the Schur complement.

Major developments were made over the years in order to build local approximations of the DtN map and accelerate the convergence of the iterative solver, including optimal [START_REF] Gander | Optimized Schwarz methods without overlap for the Helmholtz equation[END_REF], high-order [START_REF] Boubendir | A quasi-optimal non-overlapping domain decomposition algorithm for the Helmholtz equation[END_REF][START_REF] Modave | A non-overlapping domain decomposition method with high-order transmission conditions and cross-point treatment for Helmholtz problems[END_REF], Perfectly Matched Layer (PML) based [START_REF] Engquist | Sweeping preconditioner for the Helmholtz equation: moving perfectly matched layers[END_REF][START_REF] Stolk | A rapidly converging domain decomposition method for the Helmholtz equation[END_REF][START_REF] Royer | A non-overlapping domain decomposition method with perfectly matched layer transmission conditions for the Helmholtz equation[END_REF] or non-local [START_REF] Parolin | Non-overlapping domain decomposition methods with non-local transmission operators for harmonic wave propagation problems[END_REF] transmission operators. They act as a preconditioner for the global interface, Schur complement system. Let us mention that these developments are intimately related to the design of non-reflecting boundary conditions, as we will see in the article.

A global communication strategy, that allows to propagate information at a larger scale, is however necessary to achieve scalability with the number of subdomains.

Sweeping preconditionners [START_REF] Vion | Double sweep preconditioner for optimized Schwarz methods applied to the Helmholtz problem[END_REF] or coarse space corrections [START_REF] Bootland | A comparison of coarse spaces for Helmholtz problems in the high frequency regime[END_REF] could be used for that purpose. We refer to [START_REF] Gander | A class of iterative solvers for the Helmholtz equation: Factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized Schwarz methods[END_REF] for a review and general discussion of the many existing approaches. Most of these developments were applied to the Helmholtz or Maxwell equations, but less work have been devoted to problems with convection and variable coefficients with the exception of [START_REF] Tonnoir | Analysis of Schwarz methods for convected Helmholtz like equations[END_REF]. The objective of the paper is to propose a robust domain decomposition approach for industrial radiation problems governed by the Pierce equation.

The article is organized as follows. Section 2 explains how to obtain a family of high-order absorbing boundary conditions (ABCs) accounting for strong convection and variable coefficients. In particular we show robustness of the ABCs at Mach numbers close to 1 for a convex boundary shape, for which difficulties have been raised

in the literature [START_REF] Hagstrom | High-order radiation boundary conditions for the convective wave equation in exterior domains[END_REF][START_REF] Barucq | Low-order Prandtl-Glauert-Lorentz based absorbing boundary conditions for solving the convected Helmholtz equation with discontinuous Galerkin methods[END_REF]. The construction follows the pioneer approach from [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF] and extends the work initiated in [START_REF] Marchner | Construction and numerical assessment of local absorbing boundary conditions for heterogeneous time-harmonic acoustic problems[END_REF]. In Section 3, we apply the optimized Schwarz method to the Pierce operator, and propose suitable transmission operators inspired from the ABCs. We show the performance of the domain decomposition algorithm on academic problems, for which the performance is quasi-optimal for simple shaped interfaces as a direct extension to [START_REF] Boubendir | A quasi-optimal non-overlapping domain decomposition algorithm for the Helmholtz equation[END_REF]. In Section 4, we highlight the challenges of applying the domain decomposition method in more realistic settings. We address the integration of complex boundary conditions and propose a practical methodology using a second order transmission operator to handle large scale high frequency problems. We demonstrate the scalability of the algorithm on medium size problems, before showcasing its ability to solve large scale problems of industrial relevance by computing the noise radiated by the intake and the exhaust of a high by-pass ratio aircraft engine. This leads to the solution of three-dimensional finite element problems with more than a billion high-order unknowns on massively parallel, distributed memory architectures. We conclude by mentioning some perspectives and limitations.

Notations. Through all the article we use the e +ıωt time-harmonic convention, with ω the angular frequency. We introduce the vector Mach number M (x) = v 0 (x)/c 0 (x), with M = (M x , M y , M z ) T being its components in a local Cartesian coordinate system (x, y, z). The Mach number is defined to be the positive scalar quantity M = ∥M ∥, and we introduce the Prandtl-Glauert factor

β = √ 1 -M 2 .
The local wavenumber is denoted k 0 = ω/c 0 (x).

Non-reflecting boundary conditions with convection and variable

coefficients. This section is devoted to the design of non-reflecting boundary conditions for flow acoustics. We focus our study on absorbing boundary conditions, that are more direct to use in a domain decomposition context compared to other domain

This manuscript is for review purposes only.

truncation techniques such as perfectly matched layers or infinite elements.

Many time-harmonic problems are posed in an unbounded domain. It is wellknown that a Sommerfeld-like radiation condition must be satisfied at infinity in order to have a well-posed boundary value problem. Radiation conditions that take into account convection effects can be found in [START_REF] Bras | Direct and adjoint problems for sound propagation in non-uniform flows with lined and vibrating surfaces[END_REF][START_REF] Spieser | Sound propagation using an adjoint-based method[END_REF][START_REF] Bogey | Three-dimensional non-reflective boundary conditions for acoustic simulations: far field formulation and validation test cases[END_REF]]. An accurate truncation of the computational domain must be devised for volume based discretization methods, which is not obvious in the presence of spatially varying mean flow properties.

Applying a radiation condition on a finite boundary leads a non-local Dirichlet-to-Neumann (DtN) map. In this section we adopt the point of view of pseudodifferential calculus and approximate the DtN map through its symbol, that allows us to propose high-order absorbing boundary conditions. We obtain the exact DtN map for the halfspace convected Helmholtz problem, from which we deduce an approximation of the DtN map for the Pierce operator on a general convex boundary and arbitrary mean flow orientation. To implement the ABCs we use an operator localization procedure.

We conclude the section by showing numerical results of a point source radiation problem with a Mach number close to 1.

2.1. DtN map for Pierce operator. We here derive the symbol of the DtN map for the outgoing wave associated to the Pierce partial differential operator. The starting point is the formal factorization of the operator along its characteristics rays using Nirenberg's factorization theorem [START_REF] Nirenberg | Pseudodifferential Operators and Some Applications[END_REF]. The factorization guarantees that the rays are solutions to the Hamilton-Jacobi equations [START_REF] Hörmander | The Analysis of Linear Partial Differential Operators III: Pseudo-Differential Operators[END_REF], such that the characteristics carry the group velocity of the wave. The theory has been extended in [START_REF] Majda | Reflection of singularities at the boundary[END_REF] for general operators with variable coefficients. The study of the rays for the Pierce equation can be found in [START_REF] Spieser | Adjoint-based jet noise propagation model for the acoustic potential[END_REF]Chap. 2.7].

We work with the half-space domain Ω = {x ∈ R 3 ; x ≤ 0} delimited by the vertical boundary Γ = {x = 0}. We highlight that the following derivation can be extended to curved boundaries using the tangent plane approximation in a local coordinate system (n, τ 1 , τ 2 ) attached to Γ. Here n is the outward unit normal to Γ, and (τ 1 , τ 2 ) is an orthonormal basis to the tangent plane. In the half-space case we use the canonical basis of R 3 . As a preliminary step we normalize the Pierce operator

P * (ω, x, ∂ x ) = ρ 0 (1 -M 2 x ) P(ω, x, ∂ x ),
where (1 -M 2 x ) ̸ = 0 because the mean flow is assumed to be subsonic. Since the Pierce operator is a second order operator, the Nirenberg's factorization theorem for the reflected wave states that there exists two pseudodifferential operators (Λ + , Λ -) of order +1 such that

P * = ∂ x + ıΛ + ∂ x + ıΛ -+ R, (2.1)
where (Λ + , Λ -) are associated to respectively the outgoing and incoming characteristics rays, and R ∈ OPS -∞ is a regularizing pseudodifferential operator following the definition in [START_REF] Taylor | Pseudodifferential Operators[END_REF]. In order to cancel the outgoing reflected wave on the non-reflecting boundary Γ, we must satisfy the relation, for a wavefield u

∂ x u = -ıΛ + u, on Γ. (2.2)
We can then define Λ + |Γ := H 1/2 (Γ) → H -1/2 (Γ) to be the Dirichlet-to-Neumann map for the outgoing wave. We define the symbol λ + of Λ + as its partial Fourier integral

This manuscript is for review purposes only.

representation in the tangential coordinate system attached to Γ

Λ + u = R 2 λ + (τ , ξ)û(ξ)e ıτ •ξ dξ, τ ∈ Γ, (2.3)
where û(ξ) is the Fourier transform of u. We use the notation Λ + = Op(λ + ) to relate the operator to its symbol. To find the form of the operator Λ + , we expand the factorization (2.1) and identify the first and zeroth order x-derivatives terms with the Pierce operator written in the local coordinate system (x, y, z). After algebraic manipulations detailed in Appendix A, we obtain the operator equation

(1 -M 2 x ) (Λ + ) 2 + ıOp{∂ x λ + } + ı(A 1 + A 0 )(Λ + ) = (B 2 + B 1 ). (2.4)
The operators (A 1 , A 0 ) and (B 2 , B 1 ) are differential operators indexed by their degree of homogeneity with respect to (ω, ξ), because the angular frequency ω is interpreted as the symbol of the time-derivative. We have found a Ricatti-type equation for the operator Λ + , which is nothing but a reformulation of the Pierce equation as a one-way equation.

2.2.

Computation of the principal symbol. The next step is to look for the symbol λ + instead of computing Λ + directly. To do so we work with the set of classical symbols introduced by Hörmander [START_REF] Hörmander | The Analysis of Linear Partial Differential Operators III: Pseudo-Differential Operators[END_REF]. This allows to use the asymptotic

expansion λ + ∼ ∞ j=-1 λ + -j
, where each symbol is homogeneous of order (-j) with respect to (ω, ξ). The equivalence class ∼ means that the reminder is a symbol of rapid polynomial decay in ξ, associated to a smooth integral operator. Concretely we can give a notion of "convergence" for pseudodifferential operator, which here takes the form of a high frequency limit.

Moreover we have a set of calculus rules for classical symbols that allows to obtain λ + -j algebraically. Moving the operator equation (2.4) at the symbol level, we can identify the homogeneous second order terms thanks to the composition rule for pseudo-differential calculus [START_REF] Taylor | Pseudodifferential Operators[END_REF],

(1

-M 2 x )(λ + 1 ) 2 + ıσ(A 1 )λ + 1 = σ(B 2 ),
where σ( . ) denotes the symbol of the operator. Solving the polynomial of degree two for λ + gives two possibilities for the principal symbol

λ + 1 = 1 1 -M 2 x -M x (k 0 -M τ • ξ) ± k 2 0 -2k 0 M τ • ξ -β 2 |ξ| 2 , (2.5)
where M τ = (M y , M z ) T is the vector Mach number projected onto the tangential basis. The expression matches the dispersion relation for a plane wave. To select the sign in front of the square-root we can simply look that at the high frequency limit

k 0 → +∞, which gives λ + 1 → k 0 -Mx±1 1-M 2 x
. Since we target a positive group velocity for the outgoing wave, we select the positive sign and find λ + 1 → k 0 /(1 + M x ). Another point of view consists in comparing λ + 1 with the group velocity of a plane wave along the normal direction of the boundary, see [START_REF] Marchner | Non-reflecting boundary conditions and domain decomposition methods for industrial flow acoustics[END_REF]. Keeping the negative sign in front of the square-root gives the principal symbol λ - 1 associated to the incoming characteristic.

It is possible in theory to find the next terms of the symbol asymptotic expansion, as done in [START_REF] Marchner | Construction and numerical assessment of local absorbing boundary conditions for heterogeneous time-harmonic acoustic problems[END_REF] for a simplified case. In this article we simply keep the principal symbol, which amounts to truncate the expansion to the first term

λ + ≈ λ + 1 .
This manuscript is for review purposes only.

Operator representation.

The next step is to find a suitable operator related to the principal symbol λ + 1 . For the half-space problem with a uniform mean flow, it is natural to propose, using the inverse Fourier representation on Γ

Op λ + 1 = 1 1 -M 2 x -M x k 0 + ıM x (M τ • ∇ Γ ) + k 2 0 -2ık 0 (M τ • ∇ Γ ) + β 2 ∆ Γ , (2.6)
where ∇ Γ is the surface gradient and ∆ Γ is the Laplace-Beltrami operator. In this specific situation we have

Λ + = Op λ + 1 + R, R ∈ OPS -∞ .
In the general case of a non-uniform mean flow, the developments reported in Appendix A show that mean flow variations are only involved in the first order homogeneous terms for the computation of the zeroth order symbol λ + 0 . At the operator level, it implies Λ + = Op λ + 1 mod OPS 0 , in other words both operators differ by a pseudo-differential operator of order 0. The same conclusion holds when Γ is a convex curved boundary [START_REF] Antoine | Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape[END_REF].

To summarize, we have found a general high frequency approximation of the outgoing characteristic Λ + for the Pierce operator, valid for non-uniform mean flows in the tangent plane approximation. In a general curved boundary setting, we define our non-reflecting boundary condition for the outgoing wave through

∂ n u = -ıΛ + 1 u, on Γ, Λ + 1 = Op λ + 1 . (2.7)
Remark. It is worth mentioning that in some situations with variable coefficients one can exactly match the DtN symbol expansion. This occurs for example in the pure Helmholtz case where the coefficients vary only along the tangential directions.

With v 0 = (0, 0, 0) T , c 0 (x) = c 0 (y, z), ρ 0 (x) = ρ 0 (y, z), the exact DtN map for the half-space problem is

Λ + = k 2 0 + ρ 0 ∇ Γ (ρ -1 0 ∇ Γ ), on Γ.
We will show in Section 3 that a proper discretization of this operator is equivalent to a block LU factorization at the continuous level in a domain decomposition context.

2.4.

Construction of the ABCs. The operator (2.6) is still non-local. We explain in this section how to obtain local families of ABCs, following a localization procedure of the square-root operator at high frequencies.

Square-root approximation. To build a hierarchy of absorbing boundary conditions, we heavily rely on approximating the square-root function

f (z) = √ 1 + z, z ∈ C,
considering the determination along the principal branch-cut line z ∈ {-1, -∞}. We can build approximations of f thanks to a truncated Taylor or Padé expansion around z = 0, which are however not valid along the branch-cut. This is problematic since we will see that the branch-cut line is associated to evanescent modes. Following [START_REF] Milinazzo | Rational square-root approximations for parabolic equation algorithms[END_REF] we circumvent this issue by introducing a rotation of the branch-cut along the ray

R = -1 + re ı(π+α) , r > 0, α ∈ [0, -π],
such that the approximation is done on the function f α (z) = e ıα/2 e -ıα (1 + z) that we rewrite under the form

f α (z) = e ıα/2 √ 1 + z α , z α = e -ıα (1 + z) -1.
This manuscript is for review purposes only.

We can now use a classic diagonal Padé or Taylor approximation of order N ∈ N, and denote the resulting approximation respectively f PN,α and f TN,α . We compute in Figure 2.1 the absolute error

f (x) -f {T,P}N,α (x) , x ∈ R,
and show the impact of the angle α on the approximation of f . We observe that choosing α = -π/2 provides an uniform approximation along the real line except around the singular point z = -1. For a fixed order N , the choice α = 0 leads to a better approximation when z > -1 but incorrectly captures the range z < -1.

-10 -5 0 5 10 -8 Localization of the DtN map. The outgoing wave is of different nature depending on the sign of the radicand in the principal symbol (2.5), which is directly linked to microlocal regimes [START_REF] Majda | Reflection of singularities at the boundary[END_REF]. In order to capture the different regimes in a single approximation we first rewrite the operator (2.6) under the form e ıα/2 √ 1 + Z α , where Z α is a second order differential operator. We have

10 -3 10 2 x |f -f {T,P}N,α | α = 0 -10 -5 0 5 x α = -π/4 -10 -5 0 5 x α = -π/2 T1 P4 P8
Λ + 1 = -M n k 0 + ıM n (M τ • ∇ Γ ) + k 0 e ıα/2 √ 1 + Z α 1 -M 2 n , on Γ, (2.8) with Z α = e -ıα (1 + Z) -1, Z = -2ı M τ • ∇ Γ k 0 + β 2 ∆ Γ k 2 0 ,
such that Z → 0 in the high frequency limit ω → +∞. Using the square-root approximation at the operator level, we are able to approach Λ + 1 for both the evanescent and propagative parts of the spectrum when α ∈ (0, -π). This localization method allows to maintain a sparse discrete system once the operator is discretized with finite elements. The resulting local ABCs will be referred to as Taylor or Padé ABCs, called respectively ABC TN,α and ABC PN,α .

2.5. Validation for a strong uniform mean flow. The variational formulation associated to a Pierce boundary value problem is found by using Green's formula for the material derivative. For any test function v ∈ V , where V ⊂ H 1 (Ω), we obtain

Ω 1 ρ 0 ∇u • ∇v - 1 ρ 0 c 2 0 D 0 u Dt D 0 v Dt dΩ + ∂Ω v 0 • n ρ 0 c 2 0 D 0 u Dt - 1 ρ 0 ∂u ∂n v ds = 0,
This manuscript is for review purposes only.

where v is the complex conjugate, and we assumed no volume sources for conciseness.

A total flux formulation is chosen to keep the formulation symmetric, see [START_REF] Barucq | Low-order Prandtl-Glauert-Lorentz based absorbing boundary conditions for solving the convected Helmholtz equation with discontinuous Galerkin methods[END_REF] for alternative definitions. Expanding the boundary term on the non-reflecting boundary Γ reads

v 0 • n ρ 0 c 2 0 D 0 u Dt - 1 ρ 0 ∂u ∂n = ρ -1 0 ık 0 M n u + ρ -1 0 M n M τ ∂u ∂τ -ρ -1 0 1 -M 2 n ∂u ∂n . (2.9)
Replacing the normal derivative by the approximate DtN map Λ + 1 from (2.8) results

in an important simplification of the boundary integral:

Γ v 0 • n ρ 0 c 2 0 D 0 u Dt - 1 ρ 0 ∂u ∂n v ds = Γ ıe ıα/2 ρ -1 0 k 0 1 + Z α u v ds. (2.10)
We immediately note that the ABC can be implemented in the same way as in the Helmholtz case [START_REF] Boubendir | A quasi-optimal non-overlapping domain decomposition algorithm for the Helmholtz equation[END_REF][START_REF] Thierry | GetDDM: An open framework for testing optimized Schwarz methods for time-harmonic wave problems[END_REF], albeit with an adjustment of the operator Z α . We then follow the usual implementation of high-order ABCs such as presented in [START_REF] Modave | Corner treatments for high-order local absorbing boundary conditions in high-frequency acoustic scattering[END_REF]. In particular we introduce N auxiliary fields for the implementation of the Padé-based ABCs. Substituting the first order Taylor approximation Λ + 1 ≈ ABC T1,α in (2.10) leads a simple second order boundary operator

Γ ıρ -1 0 k 0 cos(α/2) + ρ -1 0 e -ıα/2 M τ ∇ Γ + ı β 2 2k 0 ∆ Γ u v ds,
which, after using Green's formula for the Laplace-Beltrami operator, can be readily implemented using the existing finite element mass and stiffness matrices. The formulation for Padé-based ABC is described in [START_REF] Marchner | Non-reflecting boundary conditions and domain decomposition methods for industrial flow acoustics[END_REF]Section 1.3.4.3]. The variational formulation is then discretized with a p-FEM conformal H 1 -basis using integrated

Legendre polynomials [START_REF] Solin | Higher-Order Finite Element Methods[END_REF].

Free field propagation. For the validation, we consider a radiating point source in a uniform mean flow located at x s , with the mean flow density and speed of sound set to unity. The exact solution is known through the Prandtl-Lorentz-Glauert transformation [START_REF] Hu | On the use of a Prandtl-Glauert-Lorentz transformation for acoustic scattering by rigid bodies with a uniform flow[END_REF] u

(x) = 1 β G (|x ′ -x ′ s |, ω ′ ) e ik 0 (M •x) β 2 , x ′ = x + (M • x) β(1 + β) M , ω ′ = ω β ,
where G is the usual Green kernel of the Helmholtz equation written in transformed coordinates. We solve the problem in two-dimensions using the finite element library

GmshFEM [START_REF] Royer | Gmsh-Fem: An Efficient Finite Element Library Based On Gmsh[END_REF]. The computational domain is a disk of unit radius, and the ABC is applied on the outer circular boundary. A Dirac source is applied at the origin. We measure the relative L 2 -error with respect to the analytical solution in the restricted domain Ω 2h = Ω\B 2h (x s ), where B 2h (x s ) is the ball of radius 2h centered at the origin of the source x s , with h the typical mesh size. The mesh is refined by a factor 5 around the point source. The polynomial order p of the shape functions is increased consistently with the Mach number in order to leverage the pollution effect. We denote d λ = 2πp k0h (1 -M ) the number of points per wavelength, that we fix to d λ = 8.

Since we deal with a purely propagative problem, there is no evanescent modes hitting the boundary Γ and we set α = 0 for the ABCs.

The real part of the solution at M = 0.95 is reported in This manuscript is for review purposes only.

1% for a Mach number close to 1 using a fourth-order Padé approximation, thanks to the combination of the high-order discretization basis and absorbing boundary condition. We have observed very little dependency to θ in our numerical results. A poor approximation of the square-root operator is only able to deal with low to moderate Mach numbers. The second order operator ABC T1,0 provides a reasonably accurate ABC for moderate Mach numbers at no supplementary numerical cost, which we keep in mind for domain decomposition. We recall that including additional curvature terms would largely improve the precision of the ABCs, especially at low Mach number [START_REF] Marchner | Non-reflecting boundary conditions and domain decomposition methods for industrial flow acoustics[END_REF]. 3. Optimized Schwarz methods for flow acoustics.

3.1. General framework. We describe in this section the optimized Schwarz formulation applied to the Pierce operator. Optimized Schwarz methods are a general framework for domain decomposition [START_REF] Gander | A class of iterative solvers for the Helmholtz equation: Factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized Schwarz methods[END_REF], and can be naturally extended to flow acoustics as initiated in [START_REF] Lieu | A nonoverlapping Schwarz domain decomposition method with high-order finite elements for flow acoustics[END_REF][START_REF] Tonnoir | Analysis of Schwarz methods for convected Helmholtz like equations[END_REF]. We hereafter use the notations and formalism following [START_REF] Antoine | Optimized Schwarz domain decomposition methods for scalar and vector Helmholtz equations[END_REF][START_REF] Boubendir | A quasi-optimal non-overlapping domain decomposition algorithm for the Helmholtz equation[END_REF]. Let us consider a non-overlapping partition of the computational domain

Ω = i Ω i , i ∈ D := {1, • • • , N dom }, with Ω i Ω j = ∅, i ̸ = j.
Each surface Γ where a boundary condition is imposed is decomposed as Γ i = Γ ∂Ω i . For a given subdomain we denote the indices of connected subdomains by D i := {j ∈ D; j ̸ = i; Σ ji ̸ = 0}. For each j ∈ D i , we define the interface that is shared between two subdomains Σ ij = Σ ji = ∂Ω i ∂Ω j . The outward oriented normal to a subdomain Ω i , respectively Ω j , is denoted n i , respectively n j . The starting point is the Pierce radiation boundary value problem in each subdomain for the unknown

u i = u| Ωi      -ρ 0 D Dt 1 ρ 2 0 c 2 0 D Dt u i + ∇ • 1 ρ0 ∇u i = 0 in Ω i , ρ -1 0 (1 -M 2 n ) (∂ ni u i + ıS i u i ) = g ij on Σ ij , ρ -1 0 (1 -M 2 n ) (∂ ni u i + ıΛ + u i ) = 0 on Γ ∞ i . (3.1)
where we have set

M 2 n := M 2 ni = M 2 nj .
The non-reflecting boundary condition on Γ ∞ i is imposed through the DtN map Λ + on Γ ∞ as described in Section 2. The operator S i is the transmission operator, which will be discussed in the next section.

The interfaces unknowns g ij are created locally for each connected subdomain. They

This manuscript is for review purposes only.

are duplicated, including at cross-points, and can be seen as Lagrange multipliers.

The interface coupling is found by eliminating the normal derivatives on a common interface Σ ij using ∂ ni u i = -∂ nj u i . We proceed as in the Helmholtz case [START_REF] Boubendir | A quasi-optimal non-overlapping domain decomposition algorithm for the Helmholtz equation[END_REF][START_REF] Thierry | GetDDM: An open framework for testing optimized Schwarz methods for time-harmonic wave problems[END_REF] and find the surface interface coupling for each subdomain i ∈ D

g ij + g ji = ıρ -1 0 (1 -M 2 n )(S i + S j )u i , on Σ ij , ∀j ∈ D i , (3.2)
where S j is the transmission operator related to Ω j . The parallel strategy is performed following an iterative procedure: each subdomain Ω i is associated to a computer process, and computes the next (n + 1) iterate as follows:

• given the interface data g

(n) ij , solve the volume problems to get u (n+1) i , • update the interface data g (n+1) ji = -g (n) ij + ıρ -1 0 (1 -M 2 n )(S i + S j )u (n+1) i .
We see that a communication pass is necessary at each iteration, in which each subdomain receives the interface data g

(n)
ij from its connected neighbourhoods and further sends g

(n+1) ji
. This strategy allows to solve volume problems of small size concurrently thanks to a direct solver, where the LU factorization is performed only once, and iterate on a naturally distributed surface problem. If the iterative method converges, such a method has the potential to handle large problems at a low memory cost.

The main concern remains to find appropriate transmission operators S i and S j to accelerate the convergence of the iterative solver for the interface problem. For that purpose we take advantage of the ABCs developed in Section 2.

Choice of transmission operator and convergence factor for the

half-space problem. The convergence factor in the two subdomains case is classically analyzed by a Fourier analysis [START_REF] Gander | Optimized Schwarz methods without overlap for the Helmholtz equation[END_REF]. We recall the analysis since it highlights important properties of our algorithm. Let us focus on the simple half-space setting with only two subdomains and a uniform mean flow with the density and speed of sound set to unity. We define the exchange operator T ji as the exchange of the outgoing trace leaving Ω i with the incoming trace coming from Ω j [START_REF] Pechstein | A unified theory of non-overlapping Robin-Schwarz methods: Continuous and discrete, including cross points[END_REF]. We write g + ij = T ji g - ji in order to highlight when the interface unknown is related to a trace with positive group velocity (+) or negative group velocity (-), which is especially important in flow acoustics. For the first exchange operator, we have g

+ ij = ∂ ni u i + ıS i u i , with
∂ ni u i = -ıΛ + u i , and g - ji = -∂ nj u j +ıS i u j , that we rewrite as g - ji = ∂ ni u i +ıS i u i . We then use ∂ ni u i = -ıΛ -u i , since (Λ + , Λ -) are respectively the outgoing and incoming DtN maps on Σ ij following Section 2. We deduce the exchange operators

T ji = (S i -Λ + )(S i -Λ -) -1 , T ij = (S j -Λ + )(S j -Λ -) -1 ,
where we have found the second exchange operator from the equality g + ji = T ij g - ij using the same reasoning. We see that the optimal transmission operator should be a good approximation of the outgoing DtN map, namely Λ + ≈ S i and Λ + ≈ S j .

However we cannot choose S j = S i as it would select the wrong mean flow orientation, and break the well-posedness of the subproblems. We choose instead S j based on the flow orientation related to the normal n j linked to Ω j , in order to guarantee that all subproblems solve the same boundary value problem. This amounts to apply the flow reversal principle. By imposing a symmetric treatment for the volume subproblems, we relocate the convection effects on the convergence factor.

Let us take a concrete example in two-dimensions and choose S i = Λ + 1 from equation (2.8). Assuming n i = (1, 0) is the outward normal related to Ω i , we have

M n = M x for S i .
Since n i = -n j , we use M n = -M x for S j by reversing the flow.

This manuscript is for review purposes only.

By substitution of the operators into the convergence factor in the partial Fourier space, we obtain

ρ(ξ) = |T ji T ij |, T ji T ij = (f N,α -f )(2M x (-k 0 + M y ξ) + f N,α -f ) (f N,α + f )(2M x (-k 0 + M y ξ) + f N,α + f ) , where f = k 2 0 -2k 0 M y ξ -β 2 ξ 2
is the symbol of the square-root operator and f N,α := f {T,P}N,α is its Taylor or Padé approximation. In the ideal case, or when f N,α is close enough to f , the convergence factor is close to zero and the algorithm converges in two iterations. In domain decomposition the convergence factor ρ acts on all the Fourier modes ξ. On the contrary to ABCs, it is here of utmost importance to select α ̸ = 0 in order to capture both the propagative and evanescent parts of the spectrum. Let us compare in Figure 3.1 the theoretical convergence factor in different situations. We show the case of a constant mean flow M = 0.7 oriented at the angles θ = {0, π/4, π/2, π} relatively to the straight interface. We refer to the transmission operators as S TN,α and S PN,α following the convention given for the ABCs in Section ), S T1,-π/2 ( ), S T1,-π/4 ( ), S P4,-π/2 ( ), S P4,-π/4 ( ).

The cases θ = {0, π} correspond to a flow oriented respectively along the positive and negative x-direction, such that M y = 0. At θ = 0 inverse upstream modes are present in the range ξ ∈ [k 0 , k 0 /β]. The square-root approximation is not valid at the two extreme values, corresponding to grazing waves, resulting in two peaks at ξ = k 0 and ξ = k 0 /β that are delimited by black lines in Figure 3.1. The convergence factor is deteriorated in this regime for all our transmission operators.

The inverse upstream regime still exists at θ = π/4 since M x > 0, but is of shorter range. When M x ≤ 0 the convergence factor is smaller than for M x > 0, for which we observe a similar behaviour to the Helmholtz case [START_REF] Boubendir | A quasi-optimal non-overlapping domain decomposition algorithm for the Helmholtz equation[END_REF]. The second order Taylor approximation with α = -π/2 provides the best trade-off between the This manuscript is for review purposes only.

attenuation of propagative and evanescent modes. It is an interesting operator since no additional cost is required for its implementation. In the Helmholtz case it shows very close performance to the second order optimized condition O02 [START_REF] Gander | Optimized Schwarz methods without overlap for the Helmholtz equation[END_REF]. For the half-space problem, Padé-based transmission operators reach arbitrary precision away from grazing modes.

3.3.

A quasi-optimal algorithm for simple partitions. The optimized Schwarz algorithm is implemented in GmshDDM, a lightweight domain decomposition library written in C++ built on top of the finite element library GmshFEM [START_REF] Royer | Gmsh-Fem: An Efficient Finite Element Library Based On Gmsh[END_REF], which is itself based on Gmsh [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in preand post-processing facilities[END_REF]. Many flow acoustics benchmarks are freely available for testing.

We here highlight the performance of the transmission operators on academic examples, and refer to Section 4 for more details regarding its usage in a HPC context.

Regarding the iterative procedure, we compute at each iteration the absolute and relative residual of the interface problem. We set our convergence criterion to a relative residual of r I = 10 -6 . When the exact solution of the problem is not available, we compute the relative L 2 -error with a usual finite element computation without domain decomposition, referred to as relative mono-domain L 2 -error. In each case we verify that the interface residual is of same order of magnitude as the mono-domain error. The problem is discretized using quadrilateral elements of typical size h and a p-FEM basis. We use d λ = 12 according to the definition provided in Section 2.5. We report in Table 3.1 the number of GMRES and Jacobi iterations to reach an interface residual of r I = 10 -6 .

Guided propagation in

Table 3.1 Number of GMRES (Jacobi) iterations at Mx = -0.7 (left) and Mx = 0.7 (right) when changing the number of subdomains. Simulation parameters: d λ = 12, p = 4, k 0 = 30 (dnc: did not converge). We observe that convergence is harder to achieve when M x > 0 due to inverse upstream modes, in particular when M x is large. The Padé conditions are very effective in this case since they can approximate the exact DtN map up to arbitrary precision on a straight interface. We observe indeed convergence in N dom iterations with a Jacobi solver. The GMRES solver accelerates the convergence of the algorithm even if some transmission conditions do not converge with the Jacobi solver.

N dom S T0,0 S T1,-π/4
This manuscript is for review purposes only.

Guided propagation without flow in a layered media. For the second example we point out that the transmission operators are also suitable for a class of problems with variable coefficients, according to the ABC performance conducted in [START_REF] Marchner | Construction and numerical assessment of local absorbing boundary conditions for heterogeneous time-harmonic acoustic problems[END_REF]. We select a normalized layered mean density and speed of sound

ρ 0 (y) = 1.25, y ∈ [0, H/2] 0.5, y ∈ [H/2, H] , c 0 (y) = 0.5, y ∈ [0, H/2] 2, y ∈ [H/2, H] .
Such a situation is encountered for turbofan exhaust acoustic radiation, see Section 4. The input is the cosine function u in (y) = cos(2πy/H), and the duct is terminated again by a PML that approximates the outgoing solution, which is validated using a mode-by-mode analysis. The key to achieve rapid convergence is to correctly localize the square-root operator. We follow an idea used in one-way modeling [START_REF] Lu | Some techniques for computing wave propagation in optical waveguides[END_REF], which consists in expressing equation ( 2.3) as follows

Λ + = ω √ 1 + Z, Z = (c -2 0 -1) + ∇ Γ (ρ -1 0 ∇ Γ )/ω 2 , (3.3)
for which we further apply the branch-cut rotation technique from Section (2.4). This operator encodes the exact DtN map for the half-space problem, and results in a correct clustering of the eigenvalues [START_REF] Marchner | Non-reflecting boundary conditions and domain decomposition methods for industrial flow acoustics[END_REF]. We report in Figure 3. With the Padé transmission operator, we observe that the number of iterations is close to the number of subdomains, which suggests that we have found the block LU factorization at the continuous level. The relative mono-domain L 2 -error indeed reaches 10 -3 % after N dom iterations, even with a Jacobi solver. The second order Taylor condition gives reasonable results. Fewer iterations are required when the second order transmission operator is used under the form (3.3).

Free field propagation. The last academic example is the point source radiation problem considered in Section 2.5. We split the domain into concentric disks of equal area, such as there are no cross-points, see Figure 3.3. We discuss the case of arbitrary partitioning in Section 4.1. We use again a PML as outgoing boundary This manuscript is for review purposes only.

condition. The mesh is refined by a factor 5 around the point source. As in Section 2.5 the shape functions polynomial order is increased steadily with the Mach number.

The mesh density is fixed to d λ = 8 for the shortest wavelength.

Conclusion.

In this section we have analyzed the convergence of the optimized Schwarz method for a class of academic problems with strong convection and variable coefficients. The proposed transmission operators perform well for unbounded waveguide and free field problems, where we have found quasi-optimal local operators in the spirit of [START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF][START_REF] Mattesi | A quasi-optimal non-overlapping domain decomposition method for two-dimensional time-harmonic elastic wave problems[END_REF]. In particular, the second order Taylor operator is interesting for its low computational cost and straightforward implementation, albeit more iterations.

4. Industrial applications and high-performance computing. We now address the development and application of the domain decomposition method in an industrial context, with the final aim to solve large scale problems on a distributed multi-core architecture. We recall that the available computer memory is the main bottleneck preventing the simulation of high frequency problems in the industry. We This manuscript is for review purposes only.

first explain how to adapt the domain decomposition methodology to industrial problems. We then briefly present two industrial benchmarks and assess the parallel scalability of our solver.

4.1.

A robust domain decomposition methodology. We have seen in the previous sections that the domain decomposition algorithm could be quasi-optimal for relatively simple partition boundaries shapes. For industrial cases we must however find a systematic way to partition a given computational domain, that might be already meshed, while ensuring a good load balancing of the subproblems unknowns.

We resort to the automatic graph partitioner METIS [START_REF] Karypis | A fast and highly quality multilevel scheme for partitioning irregular graphs[END_REF] for the generation of the subdomains, allowing a fully automatic method. It is set to minimize the size of the interface problem, and has a connectivity graph of O(N 1/d ) in d dimension, which drives the communication speed between subdomains. While the approach is very attractive at the industrial level, it also has drawbacks that affects the convergence of the iterative method. The first one is the presence of cross-points. Cross-points are a well-known theoretical and numerical issue for domain decomposition methods, and is an active topic of research [START_REF] Claeys | Robust treatment of cross-points in optimized Schwarz methods[END_REF][START_REF] Després | Corners and stable optimized domain decomposition methods for the Helmholtz problem[END_REF]. By construction of our method, we duplicate all the interface unknowns at the cross-points and do not apply any treatment on them.

In practice it has the effect to slow down the convergence of the iterative solver.

The second issue is non-smooth interface boundaries: using an automatic partitioner generates broken line interfaces, which deteriorates the quality of our transmission operators that are developed for smooth interfaces. Note that the interfaces are now cutting through the PMLs. Evanescent modes will propagate in these layers, and are naturally taken into account in our transmission conditions. The convergence to the mono-domain solution is still ensured, see also [START_REF] Modave | A non-overlapping domain decomposition method with high-order transmission conditions and cross-point treatment for Helmholtz problems[END_REF].

In order to see how such a partitioning impacts the iterative method, we reproduce in Figure 4.1 the free field propagation case from Figure 3.3 but using an arbitrary partitioning. We observe that the transmission operators still converge reasonably well, but that the Padé based conditions loose robustness with respect to the Mach number, frequency and mesh size. The second order Taylor condition here converges in a comparable number of iterations as the Padé condition. The former is hence more efficient since it has a lower computational cost, and is also easier to implement.

The choice α = -π/2 seems to be the best trade-off in such situations, where the high-order evanescent modes living on the interfaces dominate the overall error for small values of the interface residual. ), S T1,-π/2 ( ) and S P8,-π/2 ( ). In such a case the theoretical analysis is more difficult, and we view the transmis-This manuscript is for review purposes only.

sion operators as a preconditioner for the interface problem. The same observations hold in 3D, with a stronger dependency of the number of iterations with respect to the parameters k 0 , d λ and M .

To summarize, we improve the robustness of our approach at the price of complicated interface geometries, where quasi-optimal transmission conditions are much harder to design. Nevertheless, second-order conditions perform reasonably well in this context, and we will use the transmission operator S T1,-π/2 in the rest of the article.

Aeroengine intake and exhaust acoustic radiation benchmarks.

Engine noise is one of the most important acoustic source from modern commercial aircraft. On modern high by-pass ratio engines, sound is radiated by the inlet fan, by-pass duct and hot jet exhaust. These source of noise must be reduced especially during landing and take-off for health and environmental purposes.

The Pierce operator has recently been advocated to simulate such problems [START_REF] Spieser | Adjoint-based jet noise propagation model for the acoustic potential[END_REF], which has been shown to be a reasonable alternative to the Linearized Euler Equations especially in the high frequency regime [START_REF] Bras | Direct and adjoint problems for sound propagation in non-uniform flows with lined and vibrating surfaces[END_REF]. A preliminary step is to compute a steady mean flow for a given flight configuration, which is performed by an external solver.

We consider the approach and sideline configurations, representing respectively the landing and take-off flight phases, to simulate respectively the sound propagating in a hot jet exhaust and turbofan intake. For the intake problem [START_REF] Gabard | Adaptive, high-order finite-element method for convected acoustics[END_REF], a compressible, potential mean flow is solved from a non-linear Poisson problem [START_REF] Gelder | Solution of the compressible flow equations[END_REF]. The computed flow presents smooth variations but shows a strong Mach number reaching values up to M ≈ 0.8 in the throat of the nacelle, see Figure 4.2a. For the exhaust problem [START_REF] Hamiche | Hybrid numerical model for acoustic propagation through sheared flows[END_REF], a Reynolds-averaged Navier-Stokes (RANS) simulation was run.The mean flow is no more potential due to the mixing of the cold and hot jets. The Mach number shown in Figure 4.2b is moderate but the mean speed of sound and density exhibit strong variations, taking the form of a layered profile as presented in Figure 3.2. In order to avoid acoustic instabilities the mean flow boundary layer is removed as a pre-processing step [START_REF] Williamschen | Impact of the mean flow representation on DGM simulations of flow acoustics[END_REF]. Finally the pre-computed quantities c 0 , ρ 0 , v 0 are mapped on the acoustic mesh using a simple linear interpolation. modal order to be used at the bpf is given by the number of fan blades following the Tyler-Sofrin analogy [START_REF] Tyler | Axial flow compressor noise studies[END_REF]. An acoustic lining treatment is set on the walls of the nacelle duct. It models the acoustic absorption of a micro perforated plate panel backed by an air cavity, and is described by a complex-valued acoustic impedance Z(ω). The impedance is described by the Ingard-Myers boundary condition [START_REF] Myers | On the acoustic boundary condition in the presence of flow[END_REF], and implemented following [START_REF] Eversman | The boundary condition at an impedance wall in a non-uniform duct with potential mean flow[END_REF][START_REF] Bras | Direct and adjoint problems for sound propagation in non-uniform flows with lined and vibrating surfaces[END_REF]. In practice the acoustic liner has the effect to broaden the acoustic spectrum and creates back-reflected waves [START_REF] Rienstra | An Introduction to Acoustics[END_REF]. Finally a cylindrical version of the perfectly matched layer [START_REF] Marchner | Stable perfectly matched layers with Lorentz transformation for the convected Helmholtz equation[END_REF] truncates the exterior domain. A PML is also used for the input boundary condition in order to both impose the annular mode and absorb back-reflected waves.

For the exhaust problem we impose a plane wave in the heated core jet, using again a PML on the reflected field. A cylindrical version of the PML truncates the exterior domain. Both benchmarks have been validated in a two-dimensional, axi-symmetric setting, which are used for a qualitative comparison with the three-dimensional cases.

We refer to [START_REF] Gabard | Adaptive, high-order finite-element method for convected acoustics[END_REF] and [START_REF] Iob | Frequency-domain linearized euler model for turbomachinery noise radiation through engine exhaust[END_REF][START_REF] Tester | TURNEX: turbomachinery noise radiation through the engine exhaust[END_REF] for more details on the benchmarks and their applications.

4.3.

Solving large scale problems. In this last section we discuss the scalability of our solver and illustrate its capability on high frequency problems.

Parallel algorithm and scalability. Let us describe the parallel procedure used in our domain decomposition library GmshDDM. All the computations are performed in double-precision. The high-level algorithm is summarized in Algorithm 4.1. We assign one subdomain to a computer process, and use multi-threading per subdomain.

After the assembly, the LU factorization is performed once by MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF], where we use the default sparse solver. We enter the iterative solver where we rely on the PETSc GMRES implementation. We have modified the matrix-vector product thanks to a matrix free PETSc shell [START_REF] Balay | PETSc/TAO users manual[END_REF]. In our custom matrix-vector product we i) receive and send at each iteration the local interface unknowns g ij to the connected subdomains using asynchronous MPI communication, and ii) solve the subdomains volume and surface problems using the initial LU factorization. We have measured the weak scalability of our solver for a fixed number of iterations, and a problem size of roughly 800 000 unknowns per subdomain. We observed 80% scalability up to 1024 MPI processes for a total of 16 384 cores. One of the scalability limitation is the load balancing of the unknowns by the graph partitioner. With 1024 subdomains we observed a 40% relative memory increase in the peak memory usage over all the processes.

For a small number of large subdomains, the computational most intensive task is the time and memory used for the LU factorization, which scales as O(m 2 N ) with m the bandwidth of the matrix and N its size. On the contrary, with a large number of This manuscript is for review purposes only.

small subdomains the most computational intensive task is the iterative solver, since the number of iterations scales as O(N 1/3 dom ) in 3D. However the memory footprint of the iterative solver is low since it acts on a lower dimensional problem. An efficient computation relies on a trade-off between the total number of subdomains and the size of each subdomain, which is driven by the available computational resources.

A high frequency solver. We first solve a small and moderate size problem for respectively the intake and exhaust benchmarks, containing respectively around 25 and 50 wavelengths in the domain. The total problem size is respectively 10 and 86 millions of unknowns using basis functions of polynomial order p = 4, see Table 4.1. The intake case is solved at the first blade passing frequency f = 1300 Hz for the annular mode (24, 0) imposed on the fan face, and the exhaust is solved at the frequency f = 2 × 7497 Hz for the plane wave mode (0, 0) imposed at the end of the heated core duct. We solve both problems and vary the number of MPI processes.

The number of iterations and maximum peak memory usage over the MPI processes are reported in Figure 4.3. We have added estimations for the exhaust problem when computational resources were too much demanding, in particular the memory estimation was performed during the MUMPS analysis phase. A direct solver for the exhaust problem would require around 10 Tb of memory for an in-core factorization. For a low number of MPI processes, we observe a superlinear scaling in terms of computational time and memory, since the most demanding task is the LU factorization. For example the exhaust problem was solved in 7 minutes with 256 MPIs with a peak memory usage of 16.1 Gb per process, highlighting the speed and low-memory footprint of the method. Beyond a certain number of subdomains the volume subproblems are too small and there is no further gain from domain decomposition, which This manuscript is for review purposes only. can be noticed from the peak memory consumption in Figure 4.3b. The peak memory becomes bounded by the pre-processing and assembly phases, and the solving time is bounded by the iterative solver. Note that with 1024 MPIs the interface problem for the exhaust case totalizes 12 millions unknowns, but a memory footprint of less than 1 Gb and a timing of 0.4 seconds per iteration on eight CPU threads.

We continue and attempt to solve larger size problems. We set the intake problem at 5 times the blade passing frequency, f = 6500 Hz, for the input mode (48, 0) using basis functions of order p = 6, corresponding to a three-dimensional simulation of roughly 125×100×100 wavelengths in the computational domain. We split the domain into 1024 subdomains, such that each subdomain contains approximately one million unknowns. The global mesh is composed of 30 millions tetrahedra. The problem size and timing are shown in Table 4. computational domain is partitioned into 4096 subdomains, each containing approximately 320 000 unknowns. The algorithm reached a relative residual of r I = 10 -4 in 15 minutes and 555 iterations, see Table 4.3. The real part of the solution is shown in Figure 4.5. We observe guided propagation and a shortening of the wavelength in the bypass duct, as well as the wave refraction through the shear layer that impacts the acoustic radiation in the far field. As for the intake problem, we observed convergence slowdowns when increasing the input frequency, in particular for small values of the relative residual.

We ran the exhaust benchmark with 4096 MPIs for the frequencies 20, 30, 40 and 50 kHz. All reached a relative residual of r I = 10 -4 within 1000 iterations, but required respectively 1400, 3400, 6600 and 7500 GMRES iterations to reach r I = 10 -6 . Future research should investigate the introduction of a suitable coarse space to precondition the global interface problem. It could rely on interfaces plane waves [START_REF] Peng | A scalable nonoverlapping and nonconformal domain decomposition method for solving time-harmonic Maxwell equations in rˆ3[END_REF], two-level methods [START_REF] Bootland | A comparison of coarse spaces for Helmholtz problems in the high frequency regime[END_REF], or a high-frequency solver for the eikonal equation.

Conclusion and perspectives.

We have presented a hybrid direct-iterative distributed memory solver for realistic time-harmonic flow acoustic problems governed by the Pierce operator. Our solver is based on a non-overlapping optimized Schwarz method which is suitable for large scale high frequency problems that would require too much memory for direct solvers. The method is fully automatic and scales up to thousands of MPI processes, allowing to solve problems with more than a billion of high-order unknowns. These achievements were possible thanks to theoretical and This manuscript is for review purposes only. This could be coupled with advanced discretization methods such as the hybrid highorder, hybridisable discontinuous Galerkin or wave-based Treffz methods.

Finally, let us mention that domain decomposition methods are strongly linked to the physics of the problem, which makes this class of methods less robust than direct solvers. In addition domain decomposition methods require more investigations to efficiently address multiple right-hand-side.

Appendix A. Formal factorization of the Pierce operator. We present the calculations leading to the system for the operator equation (2.4). The factorization (2.1) is expanded as

P * u = ∂ 2 x + ı Λ + + Λ -∂ x + ıOp ∂ x λ + -Λ -Λ + + R, (A.1)
R ∈ OPS -∞ . We then develop the Pierce operator as

P = - 1 ρ 0 c 2 0 D 2 Dt 2 -ρ 0 v 0 • ∇ 1 ρ 2 0 c 2 0 D Dt + ∇ • 1 ρ 0 ∇ ,
and expand the double material derivative as follows

D 2 Dt 2 = -ω 2 + 2ıω v 0 • ∇ + v 0 • (v T 0 H) T + ∂V 0 ∇ , H =   ∂ 2 x ∂ 2 xy ∂ 2 xz ∂ 2 yx ∂ 2 y ∂ 2 yz ∂ 2 zx ∂ 2 zy ∂ 2 z   , ∂V 0 =   ∂ x (v x ) ∂ x (v y ) ∂ x (v z ) ∂ y (v x ) ∂ y (v y ) ∂ y (v z ) ∂ z (v x ) ∂ z (v y ) ∂ z (v z )   ,
where H is the Hessian matrix that applies to a wavefield u, and v 0 = (v x , v y , v z ) T . In that way we are able to isolate the contribution coming from the x-derivative, which is the normal derivative in a half-space setting. By identification of the first and zeroth order normal derivatives with (A.1), we form the system

ı (Λ + + Λ -) = 1 1-M 2 x (A 1 + A 0 ) iOp {∂ x λ + } -Λ -Λ + = 1 1-M 2 x (B 2 + B 1 )
,

Fig. 2 . 1 .

 21 Fig. 2.1. Absolute error between f and f {T,P}N,α along the real line.

  Fig 2.2 as well as the relative L 2 -error as a function of the Mach number. The flow angle has been fixed to θ = π/4, and we have M x = M cos(θ), M y = M sin(θ). The relative error stays below

Fig. 2 . 2 .

 22 Fig. 2.2. Numerical solution and relative ABC errors with α = 0 at k 0 = 6π as a function of the Mach number.

2. 4 .Fig. 3 . 1 .

 431 Fig. 3.1. Convergence factor for the half-space problem for different mean flow angles at M = 0.7 and k 0 = 30. Legend: S T0,0 (), S T1,-π/2 ( ), S T1,-π/4 ( ), S P4,-π/2 ( ), S P4,-π/4 ( ).

  uniform mean flow. Let us consider a multi-modal acoustic propagation in a straight waveguide geometry in a uniform mean flow only oriented along the x-direction, that is M y = 0. The waveguide is a rectangular geometry of size [0, L] × [0, H] with L = 1 and H = 0.5. It is decomposed into equally spaced subdomains separated by straight interfaces at fixed x-positions. We use a superposition of the 30 first transverse modes of unit amplitude as input boundary condition at {x = 0}. A perfectly matched layer (PML) is set as output boundary condition at {x = L}. For the moment we perform the partitioning such that the interfaces do not pass through the PML. Homogeneous Neumann boundary conditions are set on the upper and lower waveguide walls. The exact solution is available by modal superposition.

  2 the number of iterations as a function of the angular frequency to achieve convergence up to r I = 10 -6 for various transmission operators. Since we do not have an exact solution for this problem the convergence is verified by computing the relative mono-domain L 2 -error.

  (a) Real part of the solution at ω = 120. Black lines delimit the subdomains.

  Number of GMRES iterations to reach r I = 10 -6 .

Fig. 3 . 2 .

 32 Fig. 3.2. Converged numerical solution and iteration count when varying ω for the operators S T0,0 ( ), S T1,-π/4 ( ) and S P8,-π/4 ( ). Simulation parameters: N dom = 8, d λ = 12, p = 4.

4 S

 4 (a) Real part of the solution after 4 iterations at M = 0.95 using S 8,-π/2 . Black lines delimit the partitioning. P8,-π/2 (b) Number of GMRES iterations to reach r I = 10 -6 and residual history at M = 0.95 with GMRES (plain) and Jacobi (dotted) solvers.

Fig. 3 . 3 .

 33 Fig. 3.3. Numerical solution, iteration count as a function of the Mach number, and residual history for various transmission operators. Simulation parameters: k 0 = 6π, θ = π/4, d λ = 8, N dom = 5.

Figure 3 .

 3 Figure 3.3 shows the real part of the solution after the fourth GMRES iteration for the condition S P8,-π/2 with N dom = 5 at M = 0.95, θ = π/4. We immediately see the good accuracy of the numerical solution after a single forward sweep. The averaged relative L 2 -error with the exact solution is 1.8%, which is also the precision reached by the PML. We report the residual history for different transmission operators, and analyze the evolution of the number of iterations when varying the Mach number. The Padé condition shows a very fast convergence during the first iterations even with a Jacobi solver. The convergence is then slowed down due to the presence of evanescent modes on the interfaces. Using α = -π/2 instead of α = -π/4 indeed accelerates the speed of convergence. While the convergence of low order operators deteriorates at high Mach numbers when M > 0.8, Padé conditions are more robust and can almost maintain a constant number of iterations up to M = 0.95. Further tests show that Padé conditions are little sensitive to d λ and k 0 , as observed in the Helmholtz case [15].

  Variation of d λ at k 0 = 2π ( ), k 0 = 6π (), k 0 = 10π ( ) for M = 0.8 and condition S T1,-π/2 . Variation of M at k 0 = 6π, d λ = 8 with S T0,0 (

Fig. 4 . 1 .

 41 Fig. 4.1. Number of GMRES iterations to reach r I = 10 -6 for the point source free field problem with arbitrary partitioning and N dom = 16.

Fig. 4 . 2 .

 42 Fig. 4.2. Mach number in the xy-plane for the intake problem at sideline configuration and exhaust problem at approach configuration

Algorithm 4 . 1 4 .

 414 Domain decomposition algorithm for the subdomain Ω i 1. Initialization: read mesh and interpolate mean flow 2. Assembly: assemble the finite element matrix with GmshFEM, 3. Factorization: call MUMPS and run a sparse LU decomposition for the volume subproblem, Surface Assembly: assemble the surface interface problem, 5. Iterative solver: enter GMRES for the interface problem. At each iteration do until convergence: (a) receive gij and send the updated data gji to the connected subdomains, (b) compute the local matrix-vector product, 6. Post-processing

  Number of GMRES iterations to reach a relative residual of r I = 10 -6 .

Fig. 4 . 3 .

 43 Fig. 4.3. Strong scalability assessment. Iteration count and peak memory usage for the intake problem at 1300 Hz, and exhaust problem at 2×7497 Hz. Dashed lines are estimated data.

2 .

 2 The algorithm converged to a relative residual of r I = 10 -4 in 1293 iterations for a total solving time of 3 hours and 48 minutes. Note that the large input frequency increases the number of GMRES iterations because of the strong indefiniteness of the problem. The real part of the solution is shown in

Figure 4 . 4 .

 44 Figure 4.4. Strong variations of the wavelength can be observed in the throat of the nacelle.

Table 4. 2 Fig. 4 . 4 .

 244 Fig. 4.4. Real part of the solution at 5 bpf (6500 Hz), cut view of the meshed partitions, and zoom in the throat of the nacelle

Table 4. 3 Fig. 4 . 5 .

 345 Fig. 4.5. Real part of the solution at 40 kHz, cut view of the meshed partitions, and zoom in the bypass duct.

  computational contributions. On the one hand, we designed ABCs that are robust at Mach numbers close to 1 on arbitrary convex boundary shape. This led us to propose suitable transmission operators for a class of open problems with strong nonuniform mean flows. On the other hand, a high-performance framework for domain decomposition has been developed in the GmshDDM software, allowing an easy setting of industrial problems. Many improvements are already possible. At the solver level one may add p-order adaptivity, solve the interface problem on GPUs and use Block-Low Rank factorization for the subproblems. Regarding convergence, the use of overlap and PML-based transmission operator could be beneficial. A coarse space should be investigated to improve the scalability of the iterative solver. A promising perspective is to partition the computational domain with straight partitions directly through the computational mesh, and embed transmission operators as immersed boundaries.

Table 4 . 1

 41 Total problem size for the intake and exhaust problems used for the strong scalability assessment. The p-FEM order is set to 4. nnz: number of non-zeros entries.

		Frequency [Hz] unknowns	nnz	memory (est.) tetra
	Intake	1300	10M	730M	740Gb	890k
	Exhaust	2×7497	86M	6.3B	≈ 10Tb	7.7M

P. MARCHNER, H. B ÉRIOT, S. LE BRAS, X. ANTOINE AND C. GEUZAINE

This manuscript is for review purposes only.

with the differential operators

where we recall k 0 = ω/c 0 . Each operator is indexed by its degree of homogeneity with respect to (ω, ∂ y , ∂ z ). Finally we eliminate Λ -from the system and obtain the operator equation

We can readily find Λ -through