# A domain decomposition solver for large scale time-harmonic flow acoustics problems

Philippe Marchner, Hadrien Bériot, Sophie Le Bras, Xavier Antoine, Christophe Geuzaine

► To cite this version:

Philippe Marchner, Hadrien Bériot, Sophie Le Bras, Xavier Antoine, Christophe Geuzaine. A domain decomposition solver for large scale time-harmonic flow acoustics problems. SIAM Journal on Scientific Computing, 2025, 47 (2), pp.B333 - B359. 10.1137/23M1607040. hal-04254633v2

HAL Id: hal-04254633

https://hal.science/hal-04254633v2

Submitted on 13 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A DOMAIN DECOMPOSITION SOLVER FOR LARGE SCALE TIME-HARMONIC FLOW ACOUSTICS PROBLEMS *

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

**Abstract.** This article is devoted to the numerical resolution of high frequency time-harmonic flow acoustic problems. We use a substructured optimized non-overlapping Schwarz domain decomposition method as a solver in order to significantly reduce the memory footprint required by such problems. To accelerate the convergence of the iterative solver we develop suitable transmission conditions based on local approximations of the Dirichlet-to-Neumann operator, taking into account convection by strongly non-uniform mean flows. The development relies on the construction of absorbing boundary conditions through microlocal analysis and pseudo-differential calculus. We analyze the potential of the method in academic settings and subsequently propose a robust domain decomposition methodology for problems of industrial relevance modelled by the Pierce linearized acoustic operator. The algorithm is implemented in an open-source high-order finite element library. It allows to solve challenging three-dimensional problems with more than one billion high-order unknowns, by taking full advantage of modern computer architectures.

**Key words.** 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. [46, 57]), 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 [64], 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 [43]. Derived for arbitrary steady base flow, the operator is more general than the linearized potential

<footnote>*Published in SIAM Journal on Scientific Computing.</footnote>

<footnote>**Funding:** This work was performed as part of the CIFRE contract No. 2018/1845 funded by Siemens Industry Software SAS and Association Nationale de la Recherche et de la Technologie (ANRT). This research was also funded in part through the ARC grant for Concerted Research Actions (ARC WAVES 15/19-03), financed by the Wallonia-Brussels Federation of Belgium.</footnote>

<footnote>†Siemens Industry Software SAS, 107 Avenue de la République, 92320 Châtillon, France (philippe.marchner@siemens.com, sophie.le_bras@siemens.com, hadrien.beriot@siemens.com).</footnote>

<footnote>‡Université de Lorraine, CNRS, Inria, IECL, 54000 Nancy, France (xavier.antoine@univ-lorraine.fr).</footnote>

<footnote>§Université de Liège, Institut Montefiore B28, 4000 Liège, Belgium (cgeuzaine@uliege.be).</footnote>
equation from Goldstein [36] which is limited to potential homotropic mean flows. The Pierce operator has the advantage to be self-adjoint. Although its boundary conditions are in general not self-adjoint, it has been shown in [43] how to solve the adjoint problem from an existing direct solver by reversing the flow. As a result, 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 [43] for more details.

Through all this article we will consider the Pierce operator in the frequency domain. We use the $e^{+i\omega t}$ convention and define $\omega$ as the angular frequency. The Pierce operator takes the form of a convected wave operator with variable coefficients

$$ (1.1) \qquad \mathcal{P} = -\rho_0 D_{v_0} \left( \frac{1}{\rho_0^2 c_0^2} D_{v_0} \right) + \nabla \cdot \left( \frac{1}{\rho_0} \nabla \right), \quad D_{v_0} = i\omega + v_0 \cdot \nabla, $$

where the operator $D_{v_0}$ represents the material derivative along the mean flow. The quantities $\rho_0(\mathbf{x})$, $c_0(\mathbf{x})$ are respectively the spatially varying density and speed of sound of the mean flow, and $v_0(\mathbf{x})$ is the vector velocity field, which is assumed to be subsonic so that the condition $||v_0(\mathbf{x})||/c_0(\mathbf{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 [34]. 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 [7]. Moreover porous material or acoustic liners are difficult to model in the time-domain. For all these reasons an efficient and accurate frequency-domain resolution of such an operator is critical, in particular at high frequencies. In this article, we consider exterior radiation problems governed by the Pierce operator. We complete the boundary value problem by inhomogeneous and homogeneous Neumann boundary conditions to impose an acoustic velocity, Myers' boundary condition to model the impedance of acoustic liners [58, 43], and a radiation boundary condition. In this context the radiation condition can be expressed, for an outgoing radiating field $u$ [71],

$$ (1.2) \qquad \left\{ v_{g,\infty} \left( \frac{\partial}{\partial |\mathbf{x}|} + \frac{1}{|\mathbf{x}|} \right) + i\omega \right\} u(\mathbf{x}) = \mathcal{O} \left( \frac{1}{|\mathbf{x}|} \right), \quad |\mathbf{x}| \to \infty, $$

where $v_{g,\infty}$ is the wave group velocity projected along the wave propagation direction in the far-field [13], where the mean flow is assumed to be constant.

At the mathematical level, time-harmonic wave propagation problems are noto- riously difficult to solve at high frequencies [27]. The usual variational formulation of such problems are complex-valued and strongly indefinite. The strong indefinit- ness of the problem leads to inefficient iterative solvers, and makes preconditioners challenging to design [26]. As a result, direct solvers are generally advocated in the in- dustry. 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 [11]; ii) the memory requirement for direct solvers grows roughly as $\mathcal{O}(\omega^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 [69]. Order adaptivity rules have been designed [12, 10] 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 [4]. While it provides significant speed ups, memory overhead and scalability limits the applicability of direct solvers to a large number of processes [44]. 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 [3] combined with mixed-precision arithmetic [1]. 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 [30, 6] 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 convergence issue was resolved in the early 1990's by Lions [47] and Després [19] by adopting Robin transmission condition on the interfaces. An appropriate 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 [59, 25, 21], which is the continuous equivalent of the Schur complement [66]. Major developments were made over the years in order to build local approximations of DtN maps and accelerate the convergence of the iterative solver, including optimal [30], high-order [15, 56], Perfectly Matched Layer (PML) based [24, 72, 68] or non-local [18, 61] 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 [31]. Sweeping preconditioners [77] or coarse space corrections [14] could be used for that purpose. We refer to [33] 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 [32] and [45]. The objective of the paper is to propose a robust domain decomposition approach for industrial radiation problems governed by the Pierce equation. We extend our previous work [45] by providing i) a detailed analysis on how to derive and select local transmission operators for time-harmonic convected problems, and ii) a suitable numerical framework for high-performance computing.

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 [37, 9]. The construction follows the pioneer approach from [23] and extends the work initiated in [51]. 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 [15]. 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.** We introduce the vector Mach number $M(\mathbf{x}) = v_0(\mathbf{x})/c_0(\mathbf{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 $\beta = \sqrt{1 - M^2}$. The local wavenumber is denoted $k_0 = \omega/c_0(\mathbf{x})$.

**2. 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 truncation techniques such as perfectly matched layers or infinite elements.

Many time-harmonic problems are posed in an unbounded domain. It is well-known that a Sommerfeld-like radiation condition such as (1.2) must be satisfied at infinity in order to have a well-posed boundary value problem. 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.

Let us define $w$ to be the solution to (1.1) subjected to the outgoing condition (1.2) in an open exterior domain $\Omega'$. For a smooth boundary $\Gamma'$ enclosing $\Omega'$, we define a function $\varphi$ such as $\varphi = w$, on $\Gamma'$. We then define the outgoing Dirichlet-to-Neumann operator $\Lambda$ as

$$ (2.1) \qquad \Lambda\varphi = \partial_{n'}w \quad \text{on } \Gamma', $$

with $n'$ the outward unit normal to $\Gamma'$. In short, applying a radiation condition on a finite boundary leads to the definition of a Dirichlet-to-Neumann (DtN) map, which is a pseudodifferential operator of order +1 [16, 73]. 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 half-space 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. It is worth mentioning that a unified derivation of ABCs based on the Prandtl-Lorentz-Glauert has been initiated in [9]. It would be valuable, in a future work, to formulate the ABCs proposed in this article according to this framework.

**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, see Lemma 1 in [60]. The theorem has been extended in [49] for general operators with variable coefficients, which is used in the pioneer work [23]. The factorization guarantees that the rays are solutions to the Hamilton-Jacobi equations [39], such that the characteristics carry the group velocity of the wave. The study of the rays for the Pierce equation can be found in [70, Chap. 2.7].

We work with the half-space domain $\Omega = \{x \in \mathbb{R}^3; x \le 0\}$ delimited by the vertical boundary $\Gamma = \{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, \tau_1, \tau_2)$ attached to $\Gamma$. Here $n$ is the outward unit normal to $\Gamma$, and $(\tau_1, \tau_2)$ is an orthonormal basis to the tangent plane. In the half-space case we use the canonical basis of $\mathbb{R}^3$. As a preliminary step we normalize the Pierce operator

$$ \mathcal{P}^*(\omega, x, \partial_x) = \frac{\rho_0}{(1 - M_x^2)} \mathcal{P}(\omega, x, \partial_x), $$

where $(1 - M_x^2) \neq 0$ because the mean flow is assumed to be subsonic. Since the Pierce operator is a second order operator, we can apply a Nirenberg-like factorization for the reflected wave, so that there exists two pseudodifferential operators $(\Lambda^-, \Lambda^+)$ of order +1 such that

$$ (2.2) \qquad \mathcal{P}^* = (\partial_x + i\Lambda^-)(\partial_x + i\Lambda^+) + \mathcal{R}, $$

where the symbols of $(\Lambda^-, \Lambda^+)$ are associated to respectively the incoming and outgoing characteristics rays, and $\mathcal{R} \in \text{OPS}^{-\infty}$ is a regularizing pseudodifferential operator following the definition in [73]. In order to cancel the outgoing reflected wave on the non-reflecting boundary $\Gamma$, we must satisfy the relation, for a wavefield $u$

$$ (2.3) \qquad \partial_x u = -i\Lambda^+ u, \quad \text{on } \Gamma. $$

By identification with (2.1), we explicitly obtain the operator $\Lambda_{|\Gamma}^+ : H^{1/2}(\Gamma) \to H^{-1/2}(\Gamma)$ as the the Dirichlet-to-Neumann map for the outgoing wave. We define the symbol $\lambda^+$ of $\Lambda^+$ as its partial Fourier integral representation in the tangential coordinate system attached to $\Gamma$

$$ (2.4) \qquad \Lambda^+ u = \int_{\mathbb{R}^2} \lambda^+(\tau, \xi) \hat{u}(\xi) e^{-i\tau \cdot \xi} d\xi, \quad \tau \in \Gamma, $$

where $\hat{u}(\xi)$ is the Fourier transform of $u$. We use the notation $\Lambda^+ = \text{Op}(\lambda^+)$ to relate the operator to its symbol. To find the form of the operator $\Lambda^+$, we expand the factorization (2.2) 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

$$ (2.5) \qquad (1 - M_x^2) ((\Lambda^+)^2 + i \text{Op}\{\partial_x \lambda^+\}) + i(\mathcal{A}_1 + \mathcal{A}_0)(\Lambda^+) = (\mathcal{B}_2 + \mathcal{B}_1). $$
The operators ($\mathcal{A}_1, \mathcal{A}_0$) and ($\mathcal{B}_2, \mathcal{B}_1$) are differential operators indexed by their degree of homogeneity with respect to $(\omega, \xi)$, because the angular frequency $\omega$ is interpreted as the symbol of the time-derivative. We have found a Ricatti-type equation for the operator $\Lambda^+$, 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 $\lambda^+$ instead of computing $\Lambda^+$ directly. To do so we work with the set of classical symbols introduced by Hörmander [39]. This allows to use the asymptotic expansion $\lambda^+ \sim \sum_{j=-1}^{\infty} \lambda_{-j}^+$, where each symbol is homogeneous of order $(-j)$ with respect to $(\omega, \xi)$. The equivalence class $\sim$ means that the reminder is a symbol of rapid polynomial decay in $\xi$, 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 $\lambda_1^+$ algebraically. Moving the operator equation (2.5) at the symbol level, we can identify the homogeneous second order terms thanks to the composition rule for pseudo-differential calculus [73],

$$(1 - M_x^2)(\lambda_1^+) + i\sigma(\mathcal{A}_1)\lambda_1^+ = \sigma(\mathcal{B}_2),$$

where $\sigma(\cdot)$ denotes the symbol of the operator. Using Appendix A and the correspondence $\nabla_\Gamma \leftrightarrow -i\xi$ from the Fourier representation (2.4), we obtain

$$(1 - M_x^2)(\lambda_1^+) + 2M_x\lambda_1^+ (k_0 - M_\tau \cdot \xi) = (k_0 - M_\tau \cdot \xi)^2 - |\xi|^2,$$

where $M_\tau = (M_y, M_z)^T$ is the vector Mach number projected onto the tangential basis $\xi = (\xi_y, \xi_z)$. Solving the polynomial of degree two for $\lambda_1^+$ gives two possibilities for the principal symbol

$$(2.6) \quad \lambda_1^+ = \frac{1}{1 - M_x^2} \left[ -M_x(k_0 - M_\tau \cdot \xi) \pm \sqrt{(k_0 - M_\tau \cdot \xi)^2 - (1 - M_x^2)|\xi|^2} \right],$$

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 \to +\infty$, which gives $\lambda_1^+ \to k_0/(1 + M_x)$. Another point of view consists in comparing $\lambda_1^+$ with the group velocity of a plane wave along the normal direction of the boundary, see [50]. Keeping the negative sign in front of the square-root gives the principal symbol $\lambda_1^-$ associated to the incoming characteristic. It is possible in theory to find the next terms of the symbol asymptotic expansion, as done in [51] for a simplified case. In this article we simply keep the principal symbol, which amounts to truncate the expansion to the first term $\lambda^+ \approx \lambda_1^+$.

**2.3. Operator representation.** The next step is to find a suitable operator related to the principal symbol $\lambda_1^+$. For the half-space problem with a uniform mean flow, it is natural to propose, using the inverse Fourier representation on $\Gamma$

$$(2.7) \quad \mathrm{Op}(\lambda_1^+) = \frac{1}{1 - M_x^2} \left[ -M_x (k_0 - i M_\tau \cdot \nabla_\Gamma) + \sqrt{(k_0 - i M_\tau \cdot \nabla_\Gamma)^2 + (1 - M_x^2)\Delta_\Gamma} \right],$$
where $\nabla_Γ$ is the surface gradient and $\Delta_Γ$ is the Laplace-Beltrami operator. In this specific situation we have $\Lambda^+ = \text{Op}(\lambda_1^+) + \mathcal{R}$, $\mathcal{R} \in \text{OPS}^{-\infty}$. 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 $\lambda_0^+$. At the operator level, it implies $\Lambda^+ = \text{Op}(\lambda_1^+)$ mod $\text{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 [5].

To summarize, we have found a general high frequency approximation of the outgoing characteristic $\Lambda^+$ 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

$$ (2.8) \qquad \partial_n u = -i\Lambda_1^+ u, \quad \text{on } \Gamma, \quad \Lambda_1^+ = \text{Op}(\lambda_1^+). $$

**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)$, $\rho_0(x) = \rho_0(y,z)$, the exact DtN map for the half-space problem is

$$ \Lambda^+ = \sqrt{k_0^2 + \rho_0 \nabla_\Gamma (\rho_0^{-1} \nabla_\Gamma)}, \quad \text{on } \Gamma. $$

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.7) 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) = \sqrt{1+z}, \quad z \in \mathbb{C}, $$

considering the determination along the principal branch-cut line $z \in [-\infty, -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 [54] we circumvent this issue by introducing a rotation of the branch-cut along the ray $R = -1 + re^{i(\pi+\alpha)}$, $r > 0$, $\alpha \in [-\pi, 0]$, such that the approximation is done on the function $f_\alpha(z) = e^{i\alpha/2}\sqrt{e^{-i\alpha}(1+z)}$ that we rewrite under the form

$$ f_{\alpha}(z) = e^{i\alpha/2}\sqrt{1+z_{\alpha}}, \quad z_{\alpha} = e^{-i\alpha}(1+z) - 1. $$

We can now use a classic diagonal Padé or Taylor approximation of order $N \in \mathbb{N}$, and denote the resulting approximation respectively $f_{PN,\alpha}$ and $f_{TN,\alpha}$. We compute in Figure 2.1 the absolute error

$$ |f(x) - f_{\{T,P\}N,\alpha}(x)|, \quad x \in \mathbb{R}, $$

and show the impact of the angle $\alpha$ on the approximation of $f$. We observe that choosing $\alpha = -\pi/2$ provides a uniform approximation along the real line except around the singular point $z = -1$. For a fixed order $N$, the choice $\alpha = 0$ leads to a better approximation for $z > -1$ but does not capture the range $z < -1$ since $f_\alpha$ is on the branch-cut.
<figure><img src="image_2.png" /><figcaption>FIG. 2.1. Absolute error between $f$ and $f_{\{T,P\}N,\alpha}$ along the real line.</figcaption></figure>

**Localization of the DtN map.** The outgoing wave is of different nature depending on the sign of the radicand in the principal symbol (2.6), which is directly linked to microlocal regimes [49]. In order to capture the different regimes in a single approximation we first rewrite the operator (2.7) under the form $e^{i\alpha/2}\sqrt{1 + \mathcal{Z}_\alpha}$, where $\mathcal{Z}_\alpha$ is a second order differential operator. We have

$$ (2.9) \quad \Lambda_1^+ = \frac{-M_n k_0 + i M_n (M_\tau \cdot \nabla_\Gamma) + k_0 e^{i\alpha/2} \sqrt{1 + \mathcal{Z}_\alpha}}{1 - M_n^2}, \quad \text{on } \Gamma, $$

$$ \text{with } \mathcal{Z}_\alpha = e^{-i\alpha} (1 + \mathcal{Z}) - 1, \quad \mathcal{Z} = -2i \frac{M_\tau \cdot \nabla_\Gamma}{k_0} + \frac{(1 - M_n^2)\Delta_\Gamma - (M_\tau \cdot \nabla_\Gamma)^2}{k_0^2}, $$

such that $\mathcal{Z} \to 0$ in the high frequency limit $\omega \to +\infty$. Using the square-root approximation at the operator level, we are able to approach $\Lambda_1^+$ for both the evanescent and propagative parts of the spectrum when $\alpha \in (-\pi, 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,\alpha}$ and ABC$^{PN,\alpha}$.

**2.5. Validation for a strong uniform mean flow.** The variational formulation associated to a Pierce boundary value problem is derived by using Green's formula for the material derivative. For any test function $v$ in the functional space $V = \{v \in H^1(\Omega), v|_{\Gamma} \in H^1(\Gamma)\}$, we obtain

$$ \int_{\Omega} \left[ \frac{1}{\rho_0} \nabla u \cdot \overline{\nabla v} - \frac{1}{\rho_0 c_0^2} D_{v_0} u \overline{D_{v_0} v} \right] d\Omega + \int_{\partial\Omega} \left[ \frac{v_0 \cdot n}{\rho_0 c_0^2} D_{v_0} u - \frac{1}{\rho_0} \frac{\partial u}{\partial n} \right] \bar{v} ds = 0, $$

where $\bar{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 [9] for alternative definitions. Expanding the boundary term on the non-reflecting boundary $\Gamma$ reads

$$ (2.10) \quad \frac{v_0 \cdot n}{\rho_0 c_0^2} D_{v_0} u - \frac{1}{\rho_0} \frac{\partial u}{\partial n} = \rho_0^{-1} i k_0 M_n u + \rho_0^{-1} M_n M_\tau \frac{\partial u}{\partial \tau} - \rho_0^{-1} (1 - M_n^2) \frac{\partial u}{\partial n}. $$
Replacing the normal derivative by the approximate DtN map $\Lambda_1^+$ from (2.9) results in an important simplification of the boundary integral:

$$ (2.11) \qquad \int_{\Gamma} \left[ \frac{\mathbf{v}_0 \cdot \mathbf{n}}{\rho_0 c_0^2} D_{v_0} u - \frac{1}{\rho_0} \frac{\partial u}{\partial \mathbf{n}} \right] \bar{v} \, ds = \int_{\Gamma} i e^{i\alpha/2} \rho_0^{-1} k_0 \sqrt{1 + \mathcal{Z}_{\alpha} u \bar{v}} \, ds. $$

We immediately note that the ABC can be implemented in the same way as in the Helmholtz case [15, 75], albeit with an adjustment of the operator $\mathcal{Z}_{\alpha}$. We then follow the usual implementation of high-order ABCs such as presented in [55]. In particular we introduce $N$ auxiliary fields for the implementation of the Padé-based ABCs. Substituting the first order Taylor approximation $\Lambda_1^+ \approx \text{ABC}^{\text{T1},\alpha}$ in (2.11) leads to the second order boundary operator

$$ \int_{\Gamma} \left[ i \rho_0^{-1} k_0 \cos(\alpha/2) + \rho_0^{-1} e^{-i\alpha/2} \left( (M_{\tau} \cdot \nabla_{\Gamma}) + i \frac{(1-M_{n}^2)\Delta_{\Gamma} - (M_{\tau} \cdot \nabla_{\Gamma})^2}{2k_0} \right) \right] u \bar{v} \, ds, $$

which, after using Green's formula, leads to

$$ (2.12) \qquad \begin{aligned} & i \int_{\Gamma} \rho_0^{-1} k_0 \cos(\alpha/2) u \bar{v} \, ds + e^{-i\alpha/2} \int_{\Gamma} \rho_0^{-1} (M_{\tau} \cdot \nabla_{\Gamma}) u \bar{v} \, ds \\ & - i e^{-i\alpha/2} \int_{\Gamma} \frac{1-M_{n}^2}{2k_0 \rho_0} \nabla_{\Gamma} u \cdot \nabla_{\Gamma} \bar{v} \, ds + i e^{-i\alpha/2} \int_{\Gamma} \frac{1}{2k_0 \rho_0} (M_{\tau} \cdot \nabla_{\Gamma} u) (M_{\tau} \cdot \nabla_{\Gamma} \bar{v}) \, ds, \end{aligned} $$

which can be readily implemented using the existing finite element mass and stiffness matrices. The formulation for Padé-based ABC is described in [50, Section 1.3.4.3]. The variational formulation is then discretized with a p-FEM conformal $H^1$-basis using integrated Legendre polynomials [69].

**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 strong form of the problem is $\mathcal{P}u = -\delta(|x - x_s|)$, with $\delta$ the Dirac source distribution and $x_s$ the location of the source, completed with the far-field radiation condition (1.2). The exact solution is known through the Prandtl-Lorentz-Glauert transformation [40]

$$ u(x) = \frac{1}{\beta} G(|x' - x_s'|, \omega') e^{\frac{i k_0 (M \cdot x)}{\beta^2}}, \quad x' = x + \frac{(M \cdot x)}{\beta(1+\beta)} M, \quad \omega' = \frac{\omega}{\beta}, $$

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` [67]. 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 $\Omega_{2h} = \Omega \setminus 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_{\lambda} = \frac{2\pi p}{k_0 h}(1-M)$ the number of points per wavelength, that we fix to $d_{\lambda} = 8$. Since we deal with a purely propagative problem, there is no evanescent modes hitting the boundary $\Gamma$ and we set $\alpha = 0$ for the ABCs.

The real part of the solution at $M = 0.95$ is reported in 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
$\theta = \pi/4$, and we have $M_x = M \cos(\theta)$, $M_y = M \sin(\theta)$. The relative error stays below 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 $\theta$ 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. It is worth mentioning that the construction of the proposed ABC is general and can be directly applied to any boundary of convex shape. Moreover, we recall that including additional curvature terms would largely improve the precision of the ABCs, especially at low Mach number. We refer to [50] for more details.

<figure><img src="image_3.png" /><figcaption>FIG. 2.2. Numerical solution and relative ABC errors with $\alpha = 0$ at $k_0 = 6\pi$ as a function of the Mach number.</figcaption></figure>

### 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 [33], and can be naturally extended to flow acoustics as initiated in [45, 32]. We hereafter use the notations and formalism following [6, 15]. Let us consider a non-overlapping partition of the computational domain $\bar{\Omega} = \bigcup_i \bar{\Omega}_i$, $i \in D := \{1, \dots, N_{dom}\}$, with $\Omega_i \cap \Omega_j = \emptyset$, $i \neq j$. For every $i \in D$ we define the surfaces $\Gamma_i = \Gamma \cap \partial\Omega_i$ and $\Gamma_i^\infty = \Gamma^\infty \cap \partial\Omega_i$, where boundary conditions will be imposed. For $j \in D$, we define the interface that is shared between two subdomains $\bar{\Sigma}_{ij} = \bar{\Sigma}_{ji} = \partial\Omega_i \cap \partial\Omega_j$ if its interior $\Sigma_{ij}$ is not empty. For a given subdomain we denote the indices of connected subdomains by $D_i := \{j \in D; j \neq i; \Sigma_{ji} \neq 0\}$.

The outward oriented normal to a subdomain $\Omega_i$, respectively $\Omega_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|_{\Omega_i}$

$$ (3.1) \qquad \left\{ \begin{array}{ll} -\rho_0 D_{v_0} \left( \frac{1}{\rho_0^2 c_0^2} D_{v_0} u_i \right) + \nabla \cdot \left( \frac{1}{\rho_0} \nabla u_i \right) = 0 & \text{in } \Omega_i, \\ \rho_0^{-1} (1 - M_n^2) (\partial n_i u_i + i S_i u_i) = g_{ij} & \text{on } \Sigma_{ij}, \\ \rho_0^{-1} (1 - M_n^2) (\partial n_i u_i + i \Lambda^+ u_i) = 0 & \text{on } \Gamma_i^\infty. \end{array} \right. $$
where we have set $M_n^2 := M_{n_i}^2 = M_{n_j}^2$. The non-reflecting boundary condition on $\Gamma_i^\infty$ is imposed through the DtN map $\Lambda^+$ on $\Gamma_\infty$ 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 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 $\Sigma_{ij}$ using $\partial_{n_i} u_i = -\partial_{n_j} u_j$. We proceed as in the Helmholtz case [15, 75] and find the surface interface coupling for each subdomain $i \in D$

$$ (3.2) \qquad g_{ij} + g_{ji} = \iota \rho_0^{-1} (1 - M_n^2) (S_i + S_j) u_i, \quad \text{on } \Sigma_{ij}, \forall j \in D_i, $$

where $S_j$ is the transmission operator related to $\Omega_j$. The parallel strategy is performed following an iterative procedure: each subdomain $\Omega_i$ is associated to a computer process, and computes the next $(n+1)$ iterate as follows:

* given the interface data $g_{ij}^{(n)}$, solve the volume problems to get $u_i^{(n+1)}$,

* update the interface data $g_{ji}^{(n+1)} = -g_{ij}^{(n)} + \iota \rho_0^{-1} (1 - M_n^2) (S_i + S_j) u_i^{(n+1)}$.

We see that a communication pass is necessary at each iteration, in which each subdomain receives the interface data $g_{ij}^{(n)}$ from its connected neighbourhoods and further sends $g_{ji}^{(n+1)}$. 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.

## 3.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 [30]. 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 $\Omega_i$ with the incoming trace coming from $\Omega_j$ [62]. 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. Here $g_{ij}^+$ represents the amplitude of the incoming wave and $g_{ij}^-$ the amplitude of the outgoing wave. For the first exchange operator, we have $g_{ij}^+ = \partial_{n_i} u_i + \iota S_i u_i$, and select $\partial_{n_i} u_i = -\iota \Lambda^+ u_i$ to cancel the outgoing wave. For the second interface unknown we have $g_{ji}^- = \partial_{n_j} u_j + \iota S_j u_j$, that we write as $g_{ji}^- = -\iota \Lambda^- u_j + \iota S_j u_j$ in order to cancel the incoming wave. We recall that $(\Lambda^+, \Lambda^-)$ are respectively the outgoing and incoming DtN maps on $\Sigma_{ij}$ following Section 2. We deduce the exchange operators

$$ T_{ji} = (S_i - \Lambda^+)(S_j - \Lambda^-)^{-1}, \quad T_{ij} = (S_j - \Lambda^+)(S_i - \Lambda^-)^{-1}, $$

where we have found the second exchange operator from the equality $g_{ji}^+ = T_{ji}g_{ij}^-$ using the same reasoning. We see that the optimal transmission operator should be a good approximation of the outgoing DtN map, namely $\Lambda^+ \approx S_i$ and $\Lambda^+ \approx 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 $\Omega_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$ as the Taylor or Padé square-root approximation of $\Lambda_1^+$ from equation (2.9). Assuming $n_i = (1, 0)$ is the outward normal related to $\Omega_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. By denoting respectively $t_{ji}$ and $t_{ij}$ the symbols of the transmission operators $T_{ji}$ and $T_{ij}$, we can express the convergence factor in the partial Fourier space as

$$ \rho(\xi) = \sqrt{|t_{ji}(\xi)t_{ij}(\xi)|}, \quad t_{ji}t_{ij} = \frac{(f_{N,\alpha} - f)(2M_x(-k_0 + M_y\xi) + f_{N,\alpha} - f)}{(2M_x(-k_0 + M_y\xi) + f_{N,\alpha} + f)(f_{N,\alpha} + f)}, $$

where $f = \sqrt{k_0^2 - 2k_0M_y\xi - \beta^2\xi^2}$ is the symbol of the square-root operator according to equation (2.6), and $f_{N,\alpha} := f_{\{T,P\}N,\alpha}$ is its Taylor or Padé approximation. In the ideal case, or when $f_{N,\alpha}$ 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 $\rho$ acts on all the Fourier modes $\xi$. On the contrary to ABCs, it is here of utmost importance to select $\alpha \neq 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 $\theta = \{0, \pi/4, \pi/2, \pi\}$ relatively to the straight interface. We refer to the transmission operators as $S^{TN,\alpha}$ and $S^{PN,\alpha}$ following the convention given for the ABCs in Section 2.4. We show the performance of different transmission operators: the Després condition $S^{T0,0}$, the second order Taylor condition $S^{T1,\alpha}$ with $\alpha = \{-\pi/2, -\pi/4\}$ and the Padé condition $S^{P4,\alpha}$ with $N = 4$ and $\alpha = \{-\pi/2, -\pi/4\}$.

<figure><img src="image_4.png" /><figcaption>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,-\pi/2}$ (—), $S^{T1,-\pi/4}$ (---), $S^{P4,-\pi/2}$ (—), $S^{P4,-\pi/4}$ (---).</figcaption></figure>

The cases $\theta = \{0, \pi\}$ correspond to a flow oriented respectively along the positive and negative x-direction, such that $M_y = 0$. At $\theta = 0$ inverse upstream modes are
present in the range $\xi \in [k_0, k_0/\beta]$. The square-root approximation is not valid at the two extreme values, corresponding to grazing waves, resulting in two peaks at $\xi = k_0$ and $\xi = k_0/\beta$ 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 $\theta = \pi/4$ since $M_x > 0$, but is of shorter range. When $M_x \le 0$ the convergence factor is smaller than for $M_x > 0$, for which we observe a similar behaviour to the Helmholtz case [15]. The second order Taylor approximation with $\alpha = -\pi/2$ provides the best trade-off between the 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 [30]. 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` [67], which is itself based on `Gmsh` [35]. 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. We aim to achieve quasi-optimal convergence, in the sense that the required number of iterations is i) close to the depth of the connectivity graph of the subdomain partition, and ii) is quasi-independent of the spatial discretization and physical parameters such as the wavenumber or the Mach number.

**Guided propagation in 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] \times [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.

The problem is discretized using quadrilateral elements of typical size $h$ and a $p$-FEM basis. We use $d_\lambda = 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}$.

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
<caption>TABLE 3.1</caption>

Number of GMRES (Jacobi) iterations at $M_x = -0.7$ (left) and $M_x = 0.7$ (right) when changing the number of subdomains. Simulation parameters: $d_\lambda = 12, p = 4, k_0 = 30$ (dnc: did not converge).

<table><thead><tr><th>N<sub>dom</sub></th><th>S<sup>T0,0</sup></th><th>S<sup>T1,-π/4</sup></th><th>S<sup>P8,-π/4</sup></th><th>N<sub>dom</sub></th><th>S<sup>T0,0</sup></th><th>S<sup>T1,-π/4</sup></th><th>S<sup>P8,-π/4</sup></th></tr></thead><tbody><tr><td>4</td><td>28 (dnc)</td><td>19 (20)</td><td>3 (3)</td><td>4</td><td>49 (dnc)</td><td>45 (dnc)</td><td>3 (3)</td></tr><tr><td>8</td><td>61 (dnc)</td><td>43 (117)</td><td>7 (7)</td><td>8</td><td>117 (dnc)</td><td>97 (dnc)</td><td>7 (7)</td></tr><tr><td>16</td><td>155 (dnc)</td><td>87 (dnc)</td><td>15 (15)</td><td>16</td><td>268 (dnc)</td><td>197 (dnc)</td><td>17 (45)</td></tr></tbody></table>

even if some transmission conditions do not converge with the Jacobi solver.

**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 [51]. We select a normalized layered mean density and speed of sound

$$ \rho_0(y) = \begin{cases} 1.25, & y \in [0, H/2] \\ 0.5, & y \in [H/2, H] \end{cases}, \quad c_0(y) = \begin{cases} 0.5, & y \in [0, H/2] \\ 2, & y \in [H/2, H] \end{cases}. $$

Such a situation is encountered for turbofan exhaust acoustic radiation, see Section 4. The input is the cosine function $u_{in}(y) = \cos(2\pi 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 [48], which consists in expressing equation (2.3) as follows

$$ (3.3) \qquad \Lambda^{+} = \omega\sqrt{1 + \bar{Z}}, \quad \bar{Z} = (c_{0}^{-2} - 1) + \nabla_{\Gamma}(\rho_{0}^{-1}\nabla_{\Gamma})/\omega^{2}, $$

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 [50]. We report in Figure 3.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.

<figure><img src="image_5.png" /><figcaption>FIG. 3.2. Converged numerical solution and iteration count when varying $\omega$ for the operators $S^{T0,0}$ (-+–), $S^{T1,-\pi/4}$ (-*-–) and $S^{P8,-\pi/4}$ (-**–). Simulation parameters: $N_{dom} = 8, d_\lambda = 12, p = 4$.</figcaption></figure>

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 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_{\lambda} = 8$ for the shortest wavelength.

<figure><img src="image_6.png" /><figcaption>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\pi$, $\theta = \pi/4$, $d_{\lambda} = 8$, $N_{dom} = 5$.</figcaption></figure>

Figure 3.3 shows the real part of the solution after the fourth GMRES iteration for the condition $S^{P8,-\pi/2}$ with $N_{dom} = 5$ at $M = 0.95$, $\theta = \pi/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 $\alpha = -\pi/2$ instead of $\alpha = -\pi/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_{\lambda}$ and $k_0$, as observed in the Helmholtz case [15].

**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 waveg-
uide and free field problems, where we have found quasi-optimal local operators in the
spirit of [22, 53]. 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 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 [42] 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 [17, 20]. We treat cross-points by simply duplicating all the interface unknowns. 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 [56].

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 lose 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 $\alpha = -\pi/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.

In such a case the theoretical analysis is more difficult, and we view the transmission 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_\lambda$ 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^{\text{T1}, -\pi/2}$ in the rest of the
<figure><img src="image_7.png" /><figcaption>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$.</figcaption></figure>

article.

## 4.2. 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 [70], which has been shown to be a reasonable alternative to the Linearized Euler Equations especially in the high frequency regime [43]. 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 [29], a compressible, potential mean flow is solved from a non-linear Poisson problem [34]. The computed flow presents smooth variations but shows a strong Mach number reaching values up to $M \approx 0.8$ in the throat of the nacelle, see Figure 4.2a. For the exhaust problem [38], 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 [78]. Finally the pre-computed quantities $c_0, \rho_0, v_0$ are mapped on the acoustic mesh using a simple linear interpolation.

For a given mean flow, the Pierce boundary value problem can be set up. For the turbofan intake problem, tonal noise is generated at multiples of the so-called blade passing frequency (bpf). A fixed annular duct mode of azimuthal and radial orders $(m, n)$ is used to model the acoustic source on the fan face surface. The modal order to be used at the bpf is given by the number of fan blades following the Tyler-Sofrin analogy [76]. 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(\omega)$. The impedance is described by the Ingard-Myers boundary condition [58], and implemented following [28, 43]. In practice the acoustic liner has the effect to broaden the acoustic spectrum and creates back-reflected waves [65]. Finally a cylindrical version of the perfectly matched layer [52] truncates the exterior domain. A PML is
<figure><img src="image_8.png" /><figcaption>FIG. 4.2. Mach number in the xy-plane for the intake problem at sideline configuration and exhaust problem at approach configuration</figcaption></figure>

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 [29] and [41, 74] 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 [2], 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 [8]. 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.

<table><caption>**Algorithm 4.1 Domain decomposition algorithm for the subdomain $\Omega_i$**</caption><tr><td>1.</td><td><b>Initialization:</b> read mesh and interpolate mean flow</td></tr><tr><td>2.</td><td><b>Assembly:</b> assemble the finite element matrix with GmshFEM,</td></tr><tr><td>3.</td><td><b>Factorization:</b> call MUMPS and run a sparse LU decomposition for the volume subproblem,</td></tr><tr><td>4.</td><td><b>Surface Assembly:</b> assemble the surface interface problem,</td></tr><tr><td>5.</td><td><b>Iterative solver:</b> enter GMRES for the interface problem. At each iteration do until convergence:<br/>(a) receive $g_{ij}$ and send the updated data $g_{ji}$ to the connected subdomains,<br/>(b) compute the local matrix-vector product,</td></tr><tr><td>6.</td><td><b>Post-processing</b></td></tr></table>

We have measured the weak scalability of our solver for a fixed number of itera-
tions, 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^2N)$ with $m$ the bandwidth of the matrix and $N$ its size. On the contrary, with a large number of small subdomains the most computational intensive task is the iterative solver, since the number of iterations scales as $O(N_{dom}^{1/3})$ 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

<caption>TABLE 4.1</caption>

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.

<table><thead><tr><th></th><th>Frequency [Hz]</th><th>unknowns</th><th>nnz</th><th>memory (est.)</th><th>tetra</th></tr></thead><tbody><tr><td>Intake</td><td>1300</td><td>10M</td><td>730M</td><td>740Gb</td><td>890k</td></tr><tr><td>Exhaust</td><td>2×7497</td><td>86M</td><td>6.3B</td><td>≈ 10Tb</td><td>7.7M</td></tr></tbody></table>

the annular mode (24,0) imposed on the fan face, and the exhaust is solved at the frequency $f = 2 \times 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 sub-problems are too small and there is no further gain from domain decomposition, which 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
<figure><img src="image_9.png" /><figcaption>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.</figcaption></figure>

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.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. Strong variations of the wavelength can be observed in the throat of the nacelle.

<table><caption>TABLE 4.2
Total problem size and computational data for the intake problem at 5 bpf ($f = 6500$ Hz)</caption><thead><tr><th>Cores (MPI×threads)</th><th>unknowns</th><th>nnz</th><th>peak memory</th><th>pre-pro</th><th>GMRES</th><th>It</th></tr></thead><tbody><tr><td>1024×20</td><td>1.1B</td><td>167B</td><td>70Gb</td><td>24min</td><td>3h24min</td><td>1293</td></tr></tbody></table>

We also solve the exhaust benchmark set at the input frequency $f = 40$ kHz with basis function of order $p = 4$. The mesh is composed of 122 millions tetrahedra, and the shortest wavelength is discretized with approximately 6.7 degrees of freedom. The 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}$. Besides, we have observed small variations in the number of iterations depending on the METIS partitioning. Future research should investigate the introduction of a suitable coarse space to precondition the global interface problem. It could rely on interfaces plane waves [63], two-level methods [14], or a high-frequency solver for the
<figure><img src="image_10.png" /><figcaption>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</figcaption></figure>

eikonal equation.

<table><caption>TABLE 4.3
Total problem size and computational data for the exhaust problem at f = 40 kHz.</caption><thead><tr><th>Cores (MPI×threads)</th><th>unknowns</th><th>nnz</th><th>peak memory</th><th>pre-pro</th><th>GMRES</th><th>It</th></tr></thead><tbody><tr><td>4096×16</td><td>1.3B</td><td>96B</td><td>18.4Gb</td><td>1min</td><td>14min</td><td>555</td></tr></tbody></table>

**Remark.** Note that all the 3D calculations were performed with a slightly simplified version of the transmission opertor $S^{T1,-\pi/2}$, where in the last term of (2.12), the second-order cross-derivative product with $M_{\tau}$ was neglected. Numerical tests have confirmed that this simplification only led to very small differences in terms of the number of GMRES iterations, on the same order of magnitude than the differences observed when modifying the mesh partitioning.

**5. 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 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 non-uniform 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
<figure><img src="image_11.png" /><figcaption>FIG. 4.5. Real part of the solution at 40 kHz, cut view of the meshed partitions, and zoom in the bypass duct.</figcaption></figure>

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. This could be coupled with advanced discretization methods such as the hybrid high-order, 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.5). The factorization (2.2) is expanded as

$$ (A.1) \qquad \mathcal{P}^* u = \partial_x^2 + \iota (\Lambda^+ + \Lambda^-) \partial_x + \iota \mathrm{Op} \{\partial_x \lambda^+\} - \Lambda^- \Lambda^+ + \mathcal{R}, $$

$\mathcal{R} \in \mathrm{OPS}^{-\infty}$. We then develop the Pierce operator as

$$ \mathcal{P} = -\frac{1}{\rho_0 c_0^2} D_{v_0}^2 - \rho_0 \left( v_0 \cdot \nabla \left( \frac{1}{\rho_0^2 c_0^2} \right) \right) D_{v_0}^2 + \nabla \cdot \left( \frac{1}{\rho_0} \nabla \right), $$
and expand the double material derivative as follows

$$D_{v_0}^2 = -\omega^2 + 2i\omega v_0 \cdot \nabla + v_0 \cdot ((v_0^T H)^T + \partial v_0 \nabla),$$ $$H = \begin{pmatrix} \partial_x^2 & \partial_{xy}^2 & \partial_{xz}^2 \\ \partial_{yx}^2 & \partial_y^2 & \partial_{yz}^2 \\ \partial_{zx}^2 & \partial_{zy}^2 & \partial_z^2 \end{pmatrix}, \quad \partial v_0 = \begin{pmatrix} \partial_x(v_x) & \partial_x(v_y) & \partial_x(v_z) \\ \partial_y(v_x) & \partial_y(v_y) & \partial_y(v_z) \\ \partial_z(v_x) & \partial_z(v_y) & \partial_z(v_z) \end{pmatrix},$$

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

$$\left\{ \begin{array}{l} i (\Lambda^+ + \Lambda^-) = \frac{1}{1-M_x^2} (A_1 + A_0) \\ i \mathrm{Op}\{\partial_x \lambda^+\} - \Lambda^- \Lambda^+ = \frac{1}{1-M_x^2} (B_2 + B_1) \end{array} \right. ,$$

with the differential operators

$$A_1 = -2ik_0 M_x - 2M_x(M_y \partial_y + M_z \partial_z),$$

$$A_0 = -v_0 \cdot (c_0^{-2} \nabla(v_x) + \rho_0^2 \nabla(\rho_0^{-2} c_0^{-2}) v_x) + \rho_0 \partial_x (\rho_0^{-1}),$$

$$B_2 = k_0^2 - 2ik_0(M_y \partial_y + M_z \partial_z) + (1 - M_y^2) \partial_y^2 + (1 - M_z^2) \partial_z^2 - 2M_y M_z \partial_{yz}^2,$$

$$B_1 = -v_0 \cdot [c_0^{-2}(\nabla(v_y)\partial_y + \nabla(v_z)\partial_z) + \rho_0^2 \nabla(\rho_0^{-2}c_0^{-2})(v_y\partial_y + v_z\partial_z + i\omega)] \\ + \rho_0 \partial_y (\rho_0^{-1}) \partial_y + \rho_0 \partial_z (\rho_0^{-1}) \partial_z,$$

where we recall $k_0 = \omega/c_0$. Each operator is indexed by its degree of homogeneity with respect to $(\omega, \partial_y, \partial_z)$. Finally we eliminate $\Lambda^-$ from the system and obtain the operator equation

$$(1 - M_x^2)((\Lambda^+)^2 + i \mathrm{Op}\{\partial_x \lambda^+\}) + i(A_1 + A_0)(\Lambda^+) = (B_2 + B_1).$$

We can readily find $\Lambda^-$ through

$$\Lambda^- = -i \frac{1}{1 - M_x^2} (A_1 + A_0) - \Lambda^+.$$

**Reproducibility.** The source code to reproduce the results from Section 2 is available at https://gitlab.onelab.info/gmsh/fem, in the folder `examples/helmholtzFlow/freefield`. The results from Section 3 are available at https://gitlab.onelab.info/gmsh/ddm, in the folders `examples/helmholtzflow/waveguide`, `examples/helmholtz/waveguide` and `examples/helmholtzflow/freefield`. The intake problem from Section 4 is available in `examples/helmholtzflow/nacelle3D`. Instructions are provided in the corresponding `README.md` files. A Zenodo archive is also available at https://zenodo.org/records/13938470.

**Acknowledgments.** The authors would like to thank the anonymous reviewers for their critical examination of the manuscript and their insightful comments. Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11 and by the Walloon Region. The present research benefited from computational resources made available on Lucia, the Tier-1 supercomputer of the Walloon Region, infrastructure funded by the Walloon
Region under the grant agreement n°1910247. We acknowledge LUMI-BE for award-
ing this project access to the LUMI supercomputer, owned by the EuroHPC Joint
Undertaking, hosted by CSC (Finland) and the LUMI consortium through a LUMI-
BE Regular Access call. LUMI-BE is a joint effort from BELSPO (federal), SPW
Économie Emploi Recherche (Wallonia), Department of Economy, Science & Innova-
tion (Flanders) and Innoviris (Brussels). This project was provided with computer
and storage resources by GENCI at IDRIS thanks to the grant 2022-AD010613476 on
the supercomputer Jean Zay on the CSL partition. The authors would like to thank
Dr. Jérémie Gaidamour for his useful guidance and help on the supercomputer Jean
Zay.

REFERENCES

[1] P. AMESTOY, O. BOITEAU, A. BUTTARI, M. GEREST, F. JÉZÉQUEL, J.-Y. L'EXCELLENT, AND T. MARY, *Mixed precision low-rank approximations and their application to block low-rank LU factorization*, IMA Journal of Numerical Analysis, 43 (2022), pp. 2198–2227, https://doi.org/10.1093/imanum/drac037.

[2] P. AMESTOY, I. S. DUFF, J. KOSTER, AND J.-Y. L'EXCELLENT, *A fully asynchronous multifrontal solver using distributed dynamic scheduling*, SIAM Journal on Matrix Analysis and Applications, 23 (2001), pp. 15–41, https://doi.org/10.1137/S0895479899358194.

[3] P. R. AMESTOY, A. BUTTARI, J.-Y. L'EXCELLENT, AND T. MARY, *Performance and scalability of the block low-rank multifrontal factorization on multicore architectures*, ACM Trans. Math. Softw., 45 (2019), pp. 1–26, https://doi.org/10.1145/3242094.

[4] P. R. AMESTOY, I. S. DUFF, AND J.-Y. L'EXCELLENT, *Multifrontal parallel distributed symmetric and unsymmetric solvers*, Computer Methods in Applied Mechanics and Engineering, 184 (2000), pp. 501–520, https://doi.org/10.1016/S0045-7825(99)00242-X.

[5] X. ANTOINE, H. BARUCQ, AND A. BENDALI, *Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape*, Journal of Mathematical Analysis and Applications, 229 (1999), pp. 184–211, https://doi.org/10.1006/jmaa.1998.6153.

[6] X. ANTOINE AND C. GEUZAINE, *Optimized Schwarz domain decomposition methods for scalar and vector Helmholtz equations*, in *Modern Solvers for Helmholtz Problems*. Geosystems Mathematics., V. K. Lahaye D., Tang J., ed., Birkhäuser, Cham, 2017, pp. 189–213, https://doi.org/10.1007/978-3-319-28832-1_8.

[7] R. ASTLEY, R. SUGIMOTO, AND P. MUSTAFI, *Computational aero-acoustics for fan duct propagation and radiation*. Current status and application to turbofan liner optimisation, Journal of Sound and Vibration, 330 (2011), pp. 3832–3845, https://doi.org/10.1016/j.jsv.2011.03.022.

[8] S. BALAY, S. ABHYANKAR, M. ADAMS, J. BROWN, P. BRUNE, K. BUSCHELMAN, L. DALCIN, A. DENER, V. EİJKHOUT, W. GROPP, ET AL., *PETSc/TAO users manual*, 2023, https://petsc.org/release/manual/ (accessed 2023/08/21). Revision 3.19.

[9] H. BARUCQ, N. ROUXELIN, AND S. TORDEUX, *Low-order Prandtl-Glauert-Lorentz based absorbing boundary conditions for solving the convected Helmholtz equation with discontinuous Galerkin methods*, Journal of Computational Physics, 468 (2022), p. 111450, https://doi.org/10.1016/j.jcp.2022.111450.

[10] H. BÉRIOT AND G. GABARD, *Anisotropic adaptivity of the p-FEM for time-harmonic acoustic wave propagation*, Journal of Computational Physics, 378 (2019), pp. 234–256, https://doi.org/10.1016/j.jcp.2018.11.013.

[11] H. BÉRIOT, G. GABARD, AND E. PERREY-DEBAIN, *Analysis of high-order finite elements for convected wave propagation*, International Journal for Numerical Methods in Engineering, 96 (2013), pp. 665–688, https://doi.org/10.1002/nme.4559.

[12] H. BÉRIOT, A. PRINN, AND G. GABARD, *Efficient implementation of high-order finite elements for Helmholtz problems*, International Journal for Numerical Methods in Engineering, 106 (2016), pp. 213–240, https://doi.org/10.1002/nme.5172.

[13] C. BOGEY AND C. BAILLY, *Three-dimensional non-reflective boundary conditions for acoustic simulations: far field formulation and validation test cases*, Acta Acustica united with Acustica, 88 (2002), pp. 463–471.

[14] N. BOOTLAND, V. DOLEAN, P. JOLIVET, AND P.-H. TOURNIER, *A comparison of coarse spaces for Helmholtz problems in the high frequency regime*, Computers & Mathematics with
Applications, 98 (2021), pp. 239–253, https://doi.org/10.1016/j.camwa.2021.07.011.

[15] Y. BOUBENDIR, X. ANTOINE, AND C. GEUZAINE, *A quasi-optimal non-overlapping domain decomposition algorithm for the Helmholtz equation*, Journal of Computational Physics, 231 (2012), pp. 262–280, https://doi.org/10.1016/j.jcp.2011.08.007.

[16] J. CHAZARAIN AND A. PIRIOU, *Introduction to the theory of linear partial differential equations*, North-Holland, Amsterdam/New-York, 1982.

[17] X. CLAEYS AND E. PAROLIN, *Robust treatment of cross-points in optimized Schwarz methods*, Numerische Mathematik, 151 (2022), pp. 405–442, https://doi.org/10.1007/s00211-022-01288-x.

[18] F. COLLINO, S. GHANEMI, AND P. JOLY, *Domain decomposition method for harmonic wave propagation: a general presentation*, Computer methods in applied mechanics and engineering, 184 (2000), pp. 171–211, https://doi.org/10.1016/S0045-7825(99)00228-5.

[19] B. DESPRÉS, *Méthodes de Décomposition de Domaine pour les Problèmes de Propagation d'Ondes en Régime Harmonique*, Thèse, Université Paris IX Dauphine, UER Mathématiques de la Décision, (1991), http://www.theses.fr/1991PA090032.

[20] B. DESPRÉS, A. NICOLOPOULOS, AND B. THIERRY, *Corners and stable optimized domain decomposition methods for the Helmholtz problem*, Numerische Mathematik, 149 (2021), pp. 779–818, https://doi.org/10.1007/s00211-021-01251-2.

[21] V. DOLEAN, P. JOLIVET, AND F. NATAF, *An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation*, vol. 144, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2015, https://doi.org/10.1137/1.9781611974065.

[22] M. EL BOUAJAJI, B. THIERRY, X. ANTOINE, AND C. GEUZAINE, *A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations*, Journal of Computational Physics, 294 (2015), pp. 38–57, https://doi.org/10.1016/j.jcp.2015.03.041.

[23] B. ENGQUIST AND A. MAJDA, *Absorbing boundary conditions for the numerical simulation of waves*, Mathematics of Computation, 31 (1977), pp. 629–651, https://doi.org/10.2307/2005997.

[24] B. ENGQUIST AND L. YING, *Sweeping preconditioner for the Helmholtz equation: moving perfectly matched layers*, Multiscale Modeling & Simulation, 9 (2011), pp. 686–710, https://doi.org/10.1137/100804644.

[25] B. ENGQUIST AND H.-K. ZHAO, *Absorbing boundary conditions for domain decomposition*, Applied numerical mathematics, 27 (1998), pp. 341–365, https://doi.org/10.1016/S0168-9274(98)00019-1.

[26] Y. A. ERLANGGA, *Advances in iterative methods and preconditioners for the Helmholtz equation*, Archives of Computational Methods in Engineering, 15 (2008), pp. 37–66, https://doi.org/10.1007/s11831-007-9013-7.

[27] O. G. ERNST AND M. J. GANDER, *Why it is difficult to solve Helmholtz problems with classical iterative methods*, Numerical Analysis of Multiscale Problems, (2012), pp. 325–363, https://doi.org/10.1007/978-3-642-22061-6_10.

[28] W. EVERSMAN, *The boundary condition at an impedance wall in a non-uniform duct with potential mean flow*, Journal of Sound and Vibration, 246 (2001), pp. 63–69, https://doi.org/10.1006/jsvi.2000.3607.

[29] G. GABARD, H. BÉRIOT, A. PRINN, AND K. KUCUKCOSKUN, *Adaptive, high-order finite-element method for convected acoustics*, AIAA Journal, 56 (2018), pp. 3179–3191, https://doi.org/10.2514/1.J057054.

[30] M. GANDER, F. MAGOULES, AND F. NATAF, *Optimized Schwarz methods without overlap for the Helmholtz equation*, SIAM Journal on Scientific Computing, 24 (2002), pp. 38–60, https://doi.org/10.1137/S1064827501387012.

[31] M. J. GANDER AND F. KWOK, *Optimal interface conditions for an arbitrary decomposition into subdomains*, in Domain decomposition methods in science and engineering XIX, Springer, 2011, pp. 101–108, https://doi.org/10.1007/978-3-642-11304-8_9.

[32] M. J. GANDER AND A. TONNOIR, *Analysis of Schwarz Methods for Convected Helmholtz-Like Equations*, SIAM Journal on Scientific Computing, 46 (2024), pp. A1–A22, https://doi.org/10.1137/23M1560057.

[33] M. J. GANDER AND H. ZHANG, *A class of iterative solvers for the Helmholtz equation: Factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized Schwarz methods*, SIAM Review, 61 (2019), pp. 3–76, https://doi.org/10.1137/16M109781X.

[34] D. GELDER, *Solution of the compressible flow equations*, International Journal for Numerical Methods in Engineering, 3 (1971), pp. 35–43, https://doi.org/10.1002/nme.1620030107.

[35] C. GEUZAINE AND J.-F. REMACLE, *Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities*, International journal for numerical methods in engineer-
ing, 79 (2009), pp. 1309–1331, https://doi.org/10.1002/nme.2579.

[36] M. E. GOLDSTEIN, Aeroacoustics, McGraw-Hill Book Company, Inc. New York, 1976.

[37] T. HAGSTROM, S. HARIHARAN, AND D. THOMPSON, High-order radiation boundary conditions for the convective wave equation in exterior domains, SIAM Journal on Scientific Computing, 25 (2003), pp. 1088–1101, https://doi.org/10.1137/S1064827502419695.

[38] K. HAMICHE, S. LE BRAS, G. GABARD, AND H. BÉRIOT, Hybrid numerical model for acoustic propagation through sheared flows, Journal of Sound and Vibration, 463 (2019), p. 114951, https://doi.org/10.1016/j.jsv.2019.114951.

[39] L. HÖRMANDER, The Analysis of Linear Partial Differential Operators III: Pseudo-Differential Operators, Springer, 2007, https://doi.org/10.1007/978-3-540-49938-1.

[40] F. Q. HU, M. E. PIZZO, AND D. M. NARK, On the use of a Prandtl-Glauert-Lorentz transformation for acoustic scattering by rigid bodies with a uniform flow, Journal of Sound and Vibration, 443 (2019), pp. 198–211, https://doi.org/10.1016/j.jsv.2018.11.043.

[41] A. IOB, R. ARINA, AND C. SCHIPANI, Frequency-domain linearized euler model for turbomachinery noise radiation through engine exhaust, AIAA journal, 48 (2010), pp. 848–858, https://doi.org/10.2514/1.J050084.

[42] G. KARYPIS AND V. KUMAR, A fast and highly quality multilevel scheme for partitioning irregular graphs, SIAM Journal on Scientific Computing, 20 (1999), pp. 359–392, https://doi.org/10.1137/S1064827595287997.

[43] S. LE BRAS, G. GABARD, AND H. BÉRIOT, Direct and adjoint problems for sound propagation in non-uniform flows with lined and vibrating surfaces, Journal of Fluid Mechanics, 953 (2022), p. A16, https://doi.org/10.1017/jfm.2022.949.

[44] Y. LI, R. BROSSIER, AND L. MÉTIVIER, 3D frequency-domain elastic wave modeling with the spectral element method using a massively parallel direct solver, Geophysics, 85 (2020), pp. T71–T88, https://doi.org/10.1190/geo2019-0172.1.

[45] A. LIEU, P. MARCHNER, G. GABARD, H. BÉRIOT, X. ANTOINE, AND C. GEUZAINE, A non-overlapping Schwarz domain decomposition method with high-order finite elements for flow acoustics, Computer Methods in Applied Mechanics and Engineering, 369 (2020), p. 113223, https://doi.org/10.1016/j.cma.2020.113223.

[46] M. J. LIGHTHILL, On sound generated aerodynamically I. General theory, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 211 (1952), pp. 564–587, https://doi.org/10.1098/rspa.1952.0060.

[47] P.-L. LIONS, On the Schwarz alternating method. III: a variant for nonoverlapping subdomains, in Third International Symposium on Domain Decomposition Methods for Partial Differential Equations, vol. 6, SIAM Philadelphia, PA, 1990, pp. 202–223.

[48] Y. Y. LU, Some techniques for computing wave propagation in optical waveguides, Communications in Computational Physics, 1 (2006), pp. 1056–1075.

[49] A. MAJDA AND S. OSHER, Reflection of singularities at the boundary, Communications on Pure and Applied Mathematics, 28 (1975), pp. 479–499, https://doi.org/10.1002/cpa.3160280404.

[50] P. MARCHNER, Non-reflecting boundary conditions and domain decomposition methods for industrial flow acoustics, PhD thesis, Universités de Lorraine et Liège, 2022, https://www.theses.fr/2022LORR0094 (accessed 2023/08/21).

[51] P. MARCHNER, X. ANTOINE, C. GEUZAINE, AND H. BÉRIOT, Construction and numerical assessment of local absorbing boundary conditions for heterogeneous time-harmonic acoustic problems, SIAM Journal on Applied Mathematics, 82 (2022), pp. 476–501, https://doi.org/10.1137/21M1414929.

[52] P. MARCHNER, H. BÉRIOT, X. ANTOINE, AND C. GEUZAINE, Stable perfectly matched layers with Lorentz transformation for the convected Helmholtz equation, Journal of Computational Physics, 433 (2021), p. 110180, https://doi.org/10.1016/j.jcp.2021.110180.

[53] V. MATTESI, M. DARBAS, AND C. GEUZAINE, A quasi-optimal non-overlapping domain decomposition method for two-dimensional time-harmonic elastic wave problems, Journal of Computational Physics, 401 (2020), p. 109050, https://doi.org/10.1016/j.jcp.2019.109050.

[54] F. A. MILINAZZO, C. A. ZALA, AND G. H. BROOKE, Rational square-root approximations for parabolic equation algorithms, The Journal of the Acoustical Society of America, 101 (1997), pp. 760–766, https://doi.org/10.1121/1.418038.

[55] A. MODAVE, C. GEUZAINE, AND X. ANTOINE, Corner treatments for high-order local absorbing boundary conditions in high-frequency acoustic scattering, Journal of Computational Physics, 401 (2020), p. 109029, https://doi.org/10.1016/j.jcp.2019.109029.

[56] A. MODAVE, A. ROYER, X. ANTOINE, AND C. GEUZAINE, A non-overlapping domain decomposition method with high-order transmission conditions and cross-point treatment
for Helmholtz problems, Computer Methods in Applied Mechanics and Engineering, 368 (2020), p. 113162, https://doi.org/10.1016/j.cma.2020.113162.

[57] W. MOHRING, A well posed acoustic analogy based on a moving acoustic medium, Sep. 2010, https://arxiv.org/abs/1009.3766.

[58] M. MYERS, On the acoustic boundary condition in the presence of flow, Journal of Sound and Vibration, 71 (1980), pp. 429–434, https://doi.org/10.1016/0022-460X(80)90424-1.

[59] F. NATAF, F. ROGIER, AND E. DE STURLER, Optimal Interface Conditions for Domain Decomposition Methods, Internal Report 301, Ecole Polytechnique, 1994, https://hal.science/hal-02194208.

[60] L. NIRENBERG, Pseudodifferential Operators and Some Applications, in Lectures on linear partial differential equations, vol. 17, CBMS Regional Conference Series in Mathematics AMS, 1973, pp. 19–58.

[61] É. PAROLIN, Non-overlapping domain decomposition methods with non-local transmission operators for harmonic wave propagation problems, PhD thesis, Institut Polytechnique de Paris, 2020, https://www.theses.fr/2020IPPAE011 (accessed 2023/08/21).

[62] C. PECHSTEIN, A unified theory of non-overlapping Robin–Schwarz methods: Continuous and discrete, including cross points, Journal of Scientific Computing, 96 (2023), p. 60, https://doi.org/10.1007/s10915-023-02248-9.

[63] Z. PENG AND J.-F. LEE, A scalable nonoverlapping and nonconformal domain decomposition method for solving time-harmonic Maxwell equations in $\mathbb{R}^3$, SIAM Journal on Scientific Computing, 34 (2012), pp. A1266–A1295, https://doi.org/10.1137/100817978.

[64] A. D. PIERCE, Wave equation for sound in fluids with unsteady inhomogeneous flow, The Journal of the Acoustical Society of America, 87 (1990), pp. 2292–2299, https://doi.org/10.1121/1.399073.

[65] S. W. RIENSTRA AND A. HIRSCHBERG, An Introduction to Acoustics, 2004, https://www.win.tue.nl/~sjoerdr/papers/boek.pdf (accessed 2023/08/21).

[66] F.-X. ROUX, F. MAGOULES, S. SALMON, AND L. SERIES, Optimization of interface operator based on algebraic approach, in 14th International Conference on Domain Decomposition Methods, I. Herrera, D. E. Keyes, O. B. Widlund, and R. Yates, eds., National Autonomous University of Mexico (UNAM), 2002, pp. 297–304, www.ddm.org/DD14/roux.pdf.

[67] A. ROYER, E. BÉCHET, AND C. GEUZAINE, Gmsh-Fem: An Efficient Finite Element Library Based On Gmsh, in 14th World Congress on Computational Mechanics (WCCM) & EC-COMAS Congress, 2021, https://doi.org/10.23967/wccm-eccomas.2020.161.

[68] A. ROYER, C. GEUZAINE, E. BÉCHET, AND A. MODAVE, A non-overlapping domain decomposition method with perfectly matched layer transmission conditions for the Helmholtz equation, Computer Methods in Applied Mechanics and Engineering, 395 (2022), p. 115006, https://doi.org/10.1016/j.cma.2022.115006.

[69] P. SOLIN, K. SEGETH, AND I. DOLEZEL, Higher-Order Finite Element Methods, Chapman and Hall/CRC, 2003, https://doi.org/10.1201/9780203488041.

[70] E. SPIESER, Adjoint-based jet noise propagation model for the acoustic potential, PhD thesis, Université de Lyon, 2020, https://www.theses.fr/2020LYSEC043 (accessed 2023/08/21).

[71] É. SPIESER AND C. BAILLY, Sound propagation using an adjoint-based method, Journal of Fluid Mechanics, 900 (2020), p. A5, https://doi.org/10.1017/jfm.2020.469.

[72] C. STOLK, A rapidly converging domain decomposition method for the Helmholtz equation, Journal of Computational Physics, 241 (2013), pp. 240–252, https://doi.org/10.1016/j.cjp.2013.01.039.

[73] M. E. TAYLOR, Pseudodifferential Operators, Princeton University Press, 1981.

[74] B. TESTER, F. ARNOLD, S. CARO, AND S. LIDOINE, TURNEX: turbomachinery noise radiation through the engine exhaust. Publishable Final Activity Report, EU Project no, 516079 (2008).

[75] B. THIERRY, A. VION, S. TOURNIER, M. EL BOUAJAJI, D. COLIGNON, N. MARSIC, X. ANTOINE, AND C. GEUZAINE, GetDDM: An open framework for testing optimized Schwarz methods for time-harmonic wave problems, Computer Physics Communications, 203 (2016), pp. 309–330, https://doi.org/10.1016/j.cpc.2016.02.030.

[76] J. TYLER AND T. SOFRIN, Axial flow compressor noise studies, Society of Automotive Engineers Transactions, 70 (1962), pp. 309–332, https://doi.org/10.4271/620532.

[77] A. VION AND C. GEUZAINE, Double sweep preconditioner for optimized Schwarz methods applied to the Helmholtz problem, Journal of Computational Physics, 266 (2014), pp. 171–190, https://doi.org/10.1016/j.jcp.2014.02.015.

[78] M. WILLIAMSCHEN, G. GABARD, AND H. BERIOT, Impact of the mean flow representation on DGM simulations of flow acoustics, in 22nd AIAA/CEAS Aeroacoustics Conference, 2016, p. 2974, https://doi.org/10.2514/6.2016-2974.