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

HAL
open science

# A quasi-optimal coarse problem and an augmented Krylov solver for the Variational Theory of Complex Rays

L Kovalevsky, Pierre Gosselet

► To cite this version:

L Kovalevsky, Pierre Gosselet. A quasi-optimal coarse problem and an augmented Krylov solver for the Variational Theory of Complex Rays. International Journal for Numerical Methods in Engineering, 2016, 107 (11), pp.903-922. 10.1002/nme.5190 . hal-01257807

HAL Id: hal-01257807

https://hal.science/hal-01257807

Submitted on 19 Jan 2016

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A quasi-optimal coarse problem and an augmented Krylov solver for the Variational Theory of Complex Rays

L. Kovalevsky¹, P. Gosselet²

(1) Department of Engineering, University of Cambridge,
Trumpington Street, Cambridge CB2 1PZ, UK

(2) LMT-Cachan, ENS-Cachan/CNRS/Pres UniverSud Paris,
61 avenue du Président Wilson, 94235 Cachan, France

January 21, 2016

## Abstract

The Variational Theory of Complex Rays (VTCR) is an indirect Trefftz method designed to study systems governed by Helmholtz-like equations. It uses wave functions to represent the solution inside elements, which reduces the dispersion error compared to classical polynomial approaches but the resulting system is prone to be ill conditioned. This paper gives a simple and original presentation of the VTCR using the discontinuous Galerkin framework and it traces back the ill-conditioning to the accumulation of eigenvalues near zero for the formulation written in terms of wave amplitude. The core of this paper presents an efficient solving strategy that overcomes this issue. The key element is the construction of a search subspace where the condition number is controlled at the cost of a limited decrease of attainable precision. An augmented LSQR solver is then proposed to solve efficiently and accurately the complete system. The approach is successfully applied to different examples.

**Keywords:** Discoutinuous Galerkin; Helmholtz equation; Trefftz method; Variational Theory of Complex Rays; Augmented Krylov solver.

Final paper is accepted in International Journal For Numerical Methods in Engineering (doi: 10.1002/nme.5190)

## 1 Introduction

In the last decades, the use of numerical simulation techniques in design, analysis and optimization of systems has become an indispensable part of the industrial design process. The most used computer aided engineering tool is the standard Galerkin Finite Element Method (FEM [1]). It is based on the use of continuous, piecewise-polynomial shape functions supported by a mesh. It applies particularly well on the Poisson equation where a coercive formulation naturally arises.

In the case of the Helmholtz equation, the straightforward formulation is not positive, and then a fine discretization is required in order to limit the dispersion and pollution errors induced by the non-verification of the governing partial differential equations [2]. Robust FEM approaches imply to use adapted formulation and elements [3].

Alternative techniques exist based on the discontinuous Galerkin methods (DG) [4] which allow the discontinuity of the shape functions between elements, so that any type of shape function can be used. In particular, the Trefftz-DG methods use basis functions that are locally (i.e. inside each mesh element) solutions of the relevant governing partial differential equations (PDEs). In most cases, coercivity can be ensured for these methods, which moreover lead to smaller dispersion error than the finite element method (see [5] for a discussion on the error in Trefftz-DG approaches). Such approaches include, for example, the ultra-weak variational formulation (UWVF) [6], the discontinuous enrichment method [7], the wave-based method [8] and the variational theory of complex rays (VTCR) [9]. The main differences among the various Trefftz methods lie in the treatment of the boundary conditions and of the continuity conditions between elements, in the type of waves used in the admissible space and on the chosen discrete unknowns. For a given level of accuracy, all these methods lead to a substantially smaller algebraic system than one would obtain using the standard FEM. However, these methods often suffer from an ill-conditioned algebraical system (see for example [8, 10, 11, 12]) even after scaling the
different terms of the formulation [13, 7]. In order to control the condition number it was proposed in [14] to iteratively enrich the basis of wave functions until the condition number of the element matrix becomes too poor. A close idea was proposed in [12] for the VTCR.

This paper is dedicated to the study of the VTCR for acoustic problems, which is an indirect Treffz-DG method where Herglotz wave functions are used to represent the solution inside elements and where inter-elements conditions are dealt with by an anti-hermitian formulation. In Section 2 we give an original and pedagogic presentation of the VTCR and we prove that under general assumptions the formulation in terms of pressure is coercive. In Section 3 we show that the problem being set in terms of amplitudes via the (compact) Herglotz operator leads to non coercivity (even if still sign definite) because eigenvalues accumulate near zero. This phenomenon causes the bad conditioning of the discrete system. We then propose to build a search space where coercivity is restored at the price of a small decrease of the precision. In Section 4, the construction of the subspace is detailed, and it is proposed to use it as the coarse grid for an augmented Krylov solver.

The subspace that we build is quasi-optimal in the sense that, starting from an initial search space, it is a controlled approximation of the largest subspace where a chosen level of coercivity (or condition number) is preserved.

Note that this strategy is different from previous work as [14], it does not set up an upper limit for the number of waves, but selects all combinations of waves which result in pressure fields containing at least a certain amount of energy. This collection of waves spans what we call in the following the optimized subspace. The threshold used to distinguish sufficiently energetic modes from others is a parameter of the method which can be replaced by a criterion on the proportion of the total energy which should be present in the subsystem. Note that the final precision is not limited thanks to the use of the augmented LSQR solver, in that context the subspace is referred to as the coarse subspace. Two assessments are given, the first one (Section 5) possesses an analytical solution which enables us to fully understand the properties of the selected subspace and of the augmented solver; the second one (Section 6) is more realistic and enables us to illustrate the potential of the method in terms of computational performance.

# 2 The Variational Theory of Complex Rays for Helmholtz problems

## 2.1 Reference problem

Let us consider a bounded acoustic cavity $\Omega \subset \mathbb{R}^d$ filled with a fluid characterized by its speed of sound $c$, its density $\rho$ and its damping coefficient $\eta$. A source term $f$ is given in $\Omega$. The boundary $\partial\Omega$ is partitioned into three parts: $\partial_p\Omega$ where the pressure $\tilde{p}_d$ is prescribed, $\partial_v\Omega$ where the velocity $\tilde{v}_d$ is prescribed and $\partial_Z\Omega$ where a Robin condition $\tilde{h}_d$ is imposed, the impedance being written $Z$ ($Z \neq 0$ on $\partial_Z\Omega$). Assuming all excitations are time-harmonic with a given circular frequency $\omega$, the complex acoustic pressure $p$ in $\Omega$ solves the following boundary-value problem for the Helmholtz equation:

$$ \text{Find } p \text{ such that } \left. \begin{array}{lll} \Delta p + k^2 p = f & \text{in } \Omega & (a) \\ Zp + \frac{\partial p}{\partial n} = \tilde{h}_d & \text{over } \partial_Z\Omega & (b) \\ p = \tilde{p}_d & \text{over } \partial_p\Omega & (c) \\ \frac{\partial p}{\partial n} = \tilde{v}_d & \text{over } \partial_v\Omega & (d) \end{array} \right. \qquad (1) $$

where $k = (1 + i\eta)(\omega/c)$ is the wave number, $n$ is the outward normal to $\partial\Omega$, $i = \sqrt{-1}$ is the imaginary unit. For physical consideration, it is assumed that $\eta \ge 0$ and $\mathrm{Re}(Z) \ge 0$.

To apply the VTCR, an homogeneous equation is required. The reference system is then modified by the introduction of a solution $p^c$ to (1)(a) under the form $p^c(x) = \int_{\Omega} \rho(y)f(y)G(x, y) dV(y)$ where $G$ is the known Green's function<fnref n="1" />.

<footnote n="1">¹For instance, in 2D $G(x, y) = \frac{i}{4}H_0(k|x - y|)$ where $H_0$ is the zero-order Hankel function of the first kind, see [15] for general formulas.</footnote>
Setting $p^h = p - p^c$, the problem (1) can be rewritten in the following form:

$$ \text{Find } p^h \text{ such that} \quad \left\{ \begin{array}{ll} \Delta p^h + k^2 p^h = 0 & \text{in } \Omega \quad (a) \\ Zp^h + \frac{\partial p^h}{\partial n} = h_d = \tilde{h}_d - Zp^c - \frac{\partial p^c}{\partial n} & \text{over } \partial_Z \Omega \quad (b) \\ p^h = p_d = \tilde{p}_d - p^c & \text{over } \partial_p \Omega \quad (c) \\ \frac{\partial p^h}{\partial n} = v_d = \tilde{v}_d - \frac{\partial p^c}{\partial n} & \text{over } \partial_v \Omega \quad (d) \end{array} \right. \quad (2) $$

As one can see, the problem (1) with a source term, could be rewritten in an equivalent problem with modified boundary conditions and no source term. In the following we will consider a solution strategy to solve (2) and we omit superscript h

Many weak formulations of this system can be proposed, see [3] for a review, in this paper we investigate the VTCR which can be viewed as an indirect Trefftz method applied within a discontinuous-Galerkin framework.

## 2.2 One-domain VTCR formulation

Let us introduce the following space:

$$ \begin{aligned} \mathcal{V}(\Omega) &= \{u \in L^2(\Omega)/\nabla u \in L^2(\Omega)^d, \Delta u \in L^2(\Omega)\} \\ &= \{u \in H^1(\Omega)/\nabla u \in H^{\text{div}}(\Omega)\} \end{aligned} \quad (3) $$

$\mathcal{V}(\Omega)$ is an Hilbert space for the norm $\|u\|_{\mathcal{V}}^2 = \|u\|_{H^1}^2 + \|\Delta u\|_{L^2}^2$ [16]. Assuming sufficient regularity on the shape of $\Omega$, the trace and normal flux are continuous operators on $\mathcal{V}(\Omega)$ with values in $H^{\frac{1}{2}}(\partial\Omega)$ and $H^{-\frac{1}{2}}(\partial\Omega)$.

We consider the following subspace of $\mathcal{V}(\Omega)$:

$$ \mathcal{S}(\Omega) = \{p \in \mathcal{V}(\Omega) \mid \Delta p + k^2 p = 0 \text{ in } \Omega\} \quad (4) $$

$\mathcal{S}(\Omega)$ is a closed subspace of $\mathcal{V}(\Omega)$ on which the usual norm $\|u\|_{\mathcal{S}}^2 = \|\nabla u\|_{L^2}^2 + |k|^2 \|u\|_{L^2}^2$ is equivalent to the $\mathcal{V}(\Omega)$ norm and to the $H^1(\Omega)$ norm (with equivalence coefficient dependent on $k$). Moreover $\mathcal{S}(\Omega)$ is compactly embedded in $L^2(\Omega)$.

In the following we make use of the $\mathcal{S}$-norm and refer to it as the “energy norm”, though the analysis is valid for any equivalent norm (which might be more pertinent from a physical point of view).

The general one-domain VTCR formulation of the problem consists in weakly enforcing the boundary conditions as follows ($\alpha \in \mathbb{C}$ is a parameter of the formulation):

find $p \in \mathcal{S}(\Omega) / \forall q \in \mathcal{S}(\Omega)$, $a(p, q) = l(q)$ with

$$ \begin{aligned} a(p, q) &= \frac{1}{2} \int_{\partial Z \Omega} \left(Zp + \frac{\partial p}{\partial n}\right) \left(\alpha \bar{q} + \frac{\bar{\alpha}}{Z} \frac{\partial \bar{q}}{\partial n}\right) dS + \bar{\alpha} \int_{\partial_p \Omega} p \frac{\partial \bar{q}}{\partial n} dS + \alpha \int_{\partial_v \Omega} \frac{\partial p}{\partial n} \bar{q} dS \\ l(q) &= \frac{1}{2} \int_{\partial Z \Omega} h_d \left(\alpha \bar{q} + \frac{\bar{\alpha}}{Z} \frac{\partial \bar{q}}{\partial n}\right) dS + \bar{\alpha} \int_{\partial_p \Omega} p_d \frac{\partial \bar{q}}{\partial n} dS + \alpha \int_{\partial_v \Omega} v_d \bar{q} dS \end{aligned} \quad (5) $$

Clearly this formulation<fnref n="2" /> is consistent with the homogeneous reference system in the sense that the solution of (2) satisfies the weak formulation. Also $a$ is a sesquilinear form, $l$ is antilinear, both are continuous by Cauchy-Schwarz inequality and the continuity of the trace and normal flux in $\mathcal{S}(\Omega)$.

We have the following property:

$$ \begin{aligned} \operatorname{Re}(a(p, p)) &= \frac{1}{2} (\operatorname{Re}(\alpha Z)) \int_{\partial Z \Omega} \left( |p|^2 + \frac{1}{|Z|^2} \left| \frac{\partial p}{\partial n} \right|^2 \right) dS + \operatorname{Re} \left( \alpha \int_{\partial \Omega} \frac{\partial p}{\partial n} \bar{p} dS \right) \\ &= \frac{1}{2} (\operatorname{Re}(\alpha Z)) \int_{\partial Z \Omega} \left( |p|^2 + \frac{1}{|Z|^2} \left| \frac{\partial p}{\partial n} \right|^2 \right) dS \\ &\quad + \operatorname{Re}(\alpha) \|\nabla p\|_{L^2(\Omega)}^2 - \operatorname{Re}(\alpha k^2) \|p\|_{L^2(\Omega)}^2 \end{aligned} \quad (6) $$

We now discuss the well posedness of the formulation (we recall that we supposed $\operatorname{Re}(Z) \ge 0$, $\eta \ge 0$):

<footnote n="2">² To comply with the given subspace, the formulation should have been written with duality brackets in $H^{1/2}$. The integral notation does not alter computations and it is correct for the fields used in practice which are the restriction of $C^\infty(\mathbb{R}^d)$ fields, assuming sufficiently regular loads (for instance $p_d \in H^{1/2}(\partial_p \Omega)$, $v_d \in L^2(\partial_v \Omega)$, $h_d \in L^2(\partial_Z \Omega)$).</footnote>
* In the case of a negative imaginary part of the impedance and non-zero damping, the choice $\alpha = i\bar{k}$ gives a coercive formulation:

$$ \left. \begin{array}{l} \mathrm{Im}(Z) \le 0 \\ \eta > 0 \end{array} \right\} \alpha = i\bar{k} \quad \text{leads to} \quad \mathrm{Re}(a(p,p)) \ge \frac{\eta\omega}{c} \left( \|\nabla p\|_{L^2(\Omega)}^2 + |k|^2 \|p\|_{L^2(\Omega)}^2 \right) \quad (7) $$

* The domain of coercivity can be extended to impedance with small positive imaginary part if the system is sufficiently damped. Indeed, if there exists $\beta$ such that $\frac{\mathrm{Re}(Z)}{\mathrm{Im}(Z)} \ge \beta > \frac{1-\eta^2}{2\eta}$ then setting $\alpha = (1 + i\beta)$ in (6) leads to:

$$ \begin{aligned} \mathrm{Re}(a(p,p)) &= \|\nabla p\|_{L^2(\Omega)}^2 + \frac{\mathrm{Re}(Z) - \beta \mathrm{Im}(Z)}{2} \int_{\partial Z \cap \Omega} \left( |p|^2 + \frac{1}{|Z|^2} \left| \frac{\partial p}{\partial n} \right|^2 \right) dS \\ &\quad - \left( \frac{\omega^2}{c^2} \right) (1 - \eta^2 - 2\beta\eta) \|p\|_{L^2(\Omega)}^2 \\ &\ge \|\nabla p\|_{L^2(\Omega)}^2 + \left( \frac{\omega^2}{c^2} \right) \underbrace{(\eta^2 + 2\beta\eta - 1)}_{>0} \|p\|_{L^2(\Omega)}^2 \end{aligned} \quad (8) $$

The coercivity implies existence, uniqueness, and continuity with respect to the loading of the solution, and the good convergence properties of Galerkin's approximations, in particular Cea's lemma applies.

* With the minimal assumptions $\mathrm{Re}(Z) \ge 0$, $\eta \ge 0$, the simple choice $\alpha = 1$ gives a Gårding inequality:

$$ \alpha = 1 \text{ leads to } \mathrm{Re}(a(p,p)) \ge \|\nabla p\|_{L^2(\Omega)}^2 - \left(\frac{\omega}{c}\right)^2 (1 - \eta^2) \|p\|_{L^2(\Omega)}^2 \quad (9) $$

The Gårding inequality (coercivity with respect to $L^2(\Omega)$) together with the compact embedding of $S(\Omega)$ in $L^2(\Omega)$ suffices to ensure existence, uniqueness and continuity with respect to the loading when $k^2$ is not an eigenvalue of the Laplacian with the given boundary conditions; the properties of Galerkin's approximations are ensured only for sufficiently large subspaces [3].

## 2.3 Many subdomains VTCR formulation

The VTCR for many subdomains can be analysed in the framework of non-symmetric discontinuous Galerkin methods [4]. Let $\mathcal{T}$ be a partition of $\Omega$ into $N_\Omega$ non-overlapping sub-cavities $\Omega_E$ ($1 \le E \le N_\Omega$), and let $\Gamma_{E,E'} = \partial\Omega_E \cap \partial\Omega_{E'}$ be the face between subdomains $E$ and $E'$. Let $\mathcal{F}$ be the set of faces, we arbitrarily attribute an orientation to each face.

In order to treat independently the fields in each subdomain, we introduce $S^T(\Omega)$ a broken version of $S(\Omega)$:

$$ S^T(\Omega) = \{u \in L^2(\Omega) / \forall \Omega_E \in \mathcal{T}, u|_{\Omega_E} \in S(\Omega_E)\} \quad (10) $$

and we define the jump and average operators:

$$ \begin{array}{l} \forall \Gamma \in \mathcal{F}, \Gamma = \partial\Omega_E \cap \partial\Omega_{E'}, \\ \{[u]\}_\Gamma = u|_{\Omega_E} - u|_{\Omega_{E'}}, \\ \{\{u\}\}_\Gamma = \frac{u|_{\Omega_E} + u|_{\Omega_{E'}}}{2} \end{array} \quad (11) $$

We then have the following characterization of $S(\Omega)$:

$$ u \in S(\Omega) \Leftrightarrow \left\{ u \in S^T(\Omega) / \forall \Gamma \in \mathcal{F}, [[u]]_{\Gamma} = 0, \{\frac{\partial u}{\partial n}\}_{\Gamma} = 0 \right\} \quad (12) $$

In other words, fields in $S(\Omega)$ satisfy continuity and balance of normal flux conditions on the interfaces:

$$ p_E = p_{E'} \quad \text{and} \quad \frac{\partial p_E}{\partial n_E} + \frac{\partial p_{E'}}{\partial n_{E'}} = 0 \quad \text{over } \Gamma_{E,E'} $$
In the VTCR, as in discontinuous Galerkin methods, the interface conditions are introduced inside the sesquilinear form. For $(p, q) \in S^T(\Omega)^2$, we note $a_E(p, q) = a(p|_{\Omega_E}, q|_{\Omega_E})$ and $l_E(q) = l(q|_{\Omega_E})$. The VTCR for many subdomains writes:

$$ \begin{array}{l} \text{find } p \in S^T(\Omega) / \forall q \in S^T(\Omega), a_T(p, q) = l_T(q) \text{ with} \\ a_T(p, q) = \sum_E a_E(p, q) + \sum_{\Gamma \in \mathcal{F}} \int_{\Gamma} \left( \llbracket p \rrbracket_{\Gamma} \llbracket \frac{\partial \bar{q}}{\partial n} \rrbracket_{\Gamma} - \llbracket \bar{q} \rrbracket_{\Gamma} \llbracket \frac{\partial p}{\partial n} \rrbracket_{\Gamma} \right) dS \quad (13) \\ l_T(q) = \sum_E l_E(q) \end{array} $$

Thus the interface conditions are weakly imposed by an anti-hermitian formulation. This choice makes the broken formulation satisfies the same coercivity (or Gårding) inequality as the one-domain VTCR, so that no stabilization is required for the problem to be well posed. The same idea was used in Oden, Babuška and Bauman's discontinuous Galerkin formulation of the Poisson problem [17].

## 2.4 Finite-dimensional approximation space

In order to build a finite dimension approximation subspace of $S^T(\Omega)$, we use the fact that on each subcavity $\Omega_E$ the pressure can be represented by Herglotz wave functions. Let $\mathcal{C}$ be the unit sphere in $\mathbb{R}^d$ and $x_E$ a reference point located in $\Omega_E$, we define Herglotz operator $\mathcal{H}$:

$$ \begin{array}{l} \mathcal{H}_E : L^2(\mathcal{C}) \to S(\Omega_E) \\ A \mapsto \mathcal{H}_E[A] = p_A^E : x \mapsto \int_C A(s) e^{iks \cdot (x-x_E)} ds \end{array} \quad (14) $$

Under regularity assumption on the shape of $\Omega_E$ (for instance $\partial\Omega_E$ is Lipschitz and $\Omega_E$ is star-shaped with respect to the point $x_E$), range($\mathcal{H}_E$) is dense in $S(\Omega_E)$ [18, 19] so that instead of searching $p \in S^T(\Omega)$ we can seek for $(A^E) \in L^2(\mathcal{C})^{N_{\Omega}}$. Keeping terminology adopted in previous work on the VTCR, the density $A^E$ is called the amplitude distribution of $p_A^E$ since somehow the pressure is represented by a superposition of plane waves $e^{iks \cdot (x-x_E)}$ in direction $s$ with magnitude $A^E(s)$. The formulation in terms of unknown amplitudes makes the VTCR an indirect Trefftz method.

The retained strategy is then to discretize the space of amplitude distribution $L^2(\mathcal{C})^{N_{\Omega}}$. For a given subdomain $E$, the approximation subspace $\mathcal{A}^E$ of $L^2(\mathcal{C})$ of dimension $N_E$ is defined by a basis $\mathbf{A}^E = (..., A_n^E, ...)$.

The resulting subspace of pressure is denoted $S^{N_E}(\Omega^E)$, it is spanned by $\mathbf{P}_A^E = (..., \mathcal{H}_E(A_n^E), ...)$.

Various discretizations have been tested in previous implementations of the VTCR [11, 20]. To keep expressions simple, we recall them in the 2D case where $s$ only depends on one angle $\theta$: $s_{\theta} = \begin{pmatrix} \cos(\theta) \\ \sin(\theta) \end{pmatrix}$.

$D_B$: $\mathbf{A}^E$ is made out of piecewise-constant functions on $\mathcal{C}$, $\mathbf{P}_A^E$ is then a collection of wave Band functions:

$$ \mathbf{P}_A^E = \left( \dots, \int_{\frac{2\pi n}{N_E}}^{\frac{2\pi(n+1)}{N_E}} e^{iks_\theta \cdot (x-x_E)} d\theta, \dots \right) \quad (15) $$

$D_F$: $\mathbf{A}^E$ is taken as a truncated Fourier series, $\mathbf{P}_A^E$ then being a collection of Fourier wave functions ($N_E$ is assumed to be odd):

$$ \mathbf{P}_A^E = \left( \dots, \int_{-\pi}^{\pi} e^{iks_\theta} e^{iks_\theta \cdot (x-x_E)} d\theta, \dots \right) \quad (16) $$

The main advantages of this discretization is that the functions of $\mathbf{P}_A^E$ can be evaluated analytically at any point $x$. Moreover these functions form a hierarchical basis, which makes it easy to increase the degree of the approximation, in particular a control of the condition number, similar to [14], was proposed in [20].

$D_{\delta}$: One last possibility is to extend the definition of Herglotz to a larger space containing Dirac distributions$^3$. The typical choice is $\mathbf{A}^E$ constituted by $N_E$ Dirac distributions supported at angular locations $\theta_n = \frac{2\pi n}{N_E}$

<footnote>$^3$For instance $H^{-s}(\mathcal{C})$ with $s > \frac{d-1}{2}$ where $d$ is the dimension of the physical space</footnote>
(in some publications these functions are also called “rays”). $\mathbf{P}_A^E$ is then the collection of a finite number of plane waves:

$$ \mathbf{P}_A^E = (\dots, e^{i k s_{\theta_n} \cdot (x-x_E)}, \dots)_{0 \le n < N_E} \qquad (17) $$

This discretization presents the main advantage to allow analytical integration of the components of the weak formulation (13) on straight lines, it also corresponds exactly to the enrichment functions used in the Discontinuous Enrichment Method [7]. Unfortunately it does not fit exactly the theory described in the following, which is why we propose an informal extension to that discretization in section 3.5.

The above discretizations, which have their pros and cons, lead to very similar asymptotic accuracy, as seen in comparisons presented in [11, 20]. From these studies, a priori criteria have emerged to choose a sufficiently fine discretization, represented by the dimension of the search space $N_E$. In the Fourier case (16), $N_E$ can be set a priori using an energetic criterion [20]. In Dirac (17) and Band (15) cases, $N_E$ can be set a priori using a geometrical heuristic criterion

$$ N_E = \text{round}(\mu\pi R_E/\lambda) \qquad (18) $$

where $R_E$ is the characteristic diameter of $\Omega_E$, $\lambda$ is the wavelength and $\mu$ a positive real number close to 1 (see [21] eqn. 43, or [22] eqn. 3.37, and the example of section 5).

## 2.5 Discrete system

Introducing any of the previous basis in the weak formulation (13) leads to a linear system with the following block structure:

$$ \mathbf{K}\alpha = \mathbf{f}, \text{ with } \mathbf{K} = \begin{pmatrix} \mathbf{K}_{11} & \mathbf{K}_{12} & \dots & \mathbf{K}_{1N_{\Omega}} \\ -\mathbf{K}_{12}^H & \mathbf{K}_{22} & & \mathbf{K}_{2N_{\Omega}} & \dots & \mathbf{K}_{N_{\Omega}N_{\Omega}} \\ \vdots & \vdots & & \vdots & & \vdots \\ -\mathbf{K}_{1N_{\Omega}}^H & \mathbf{K}_{1N_{\Omega}} & \dots & \mathbf{K}_{1N_{\Omega}} & \dots & \mathbf{K}_{N_{\Omega}N_{\Omega}} \end{pmatrix}, \alpha = \begin{pmatrix} \alpha_1 \\ \vdots \\ \alpha_{N_{\Omega}} \end{pmatrix}, \mathbf{f} = \begin{pmatrix} \mathbf{f}_1 \\ \vdots \\ \mathbf{f}_{N_{\Omega}} \end{pmatrix} \qquad (19) $$

where exponent $H$ stands for the conjugation-transposition, and $\alpha_E$ corresponds to the vector of unknown amplitude of the chosen basis functions in the subdomain $\Omega_E$. The off-diagonal blocks correspond to the coupling between subdomains, they are zero for non-neighbors subdomains so that the system has a sparse-by-block structure. In the case where damping ensures coercivity of the formulation, the diagonal blocks are positive.

Each degree of freedom is associated with one wave function in one subdomain and it necessarily contributes to the off-diagonal coupling term. This implies that contrary to what is commonly encountered in DG methods [23, 24] there is no “internal” degree of freedom and condensation can not apply.

# 3 The Optimized space of approximation

One major problem with searching for amplitude components is that operator $\mathcal{H}_E$ is compact in $L^2(C)$: there is an accumulation of eigenvalues near zero. In other words there exist amplitudes of unit norm in $L^2(C)$ capable to create arbitrary small pressure fields in $S(\Omega_E)$:

$$ \forall \varepsilon > 0, \exists A \in L^2(C) \text{ with } \|A\|_{L^2(C)} = 1 \text{ and } \|\mathcal{H}_E[A]\|_{S(\Omega_E)} < \varepsilon \qquad (20) $$

To illustrate the consequences of this problem, let us consider the favorable case of a coercive formulation in $S^T(\Omega)$ where there exists a positive constant $Q$ such that:

$$ \forall p \in S^T(\Omega), a_T(p, p) \ge Q \|p\|_{S(\Omega)}^2 \qquad (21) $$

If we now consider the formulation in terms of amplitudes, all we can say is:

$$ \forall (A^E) \in L^2(C)^{N_{\Omega}}, a_T((\mathcal{H}_E[A^E]), (\mathcal{H}_E[A^E])) \ge 0 \qquad (22) $$

the sesquilinear form is positive but it can not be bounded from below by the norm of the amplitudes $\sum_E \|A^E\|_{L^2(C)}^2$, the formulation in amplitude is not coercive. Moreover, the discrete system is likely to be poorly conditioned: as soon as the discretization space is large enough, it is probable that it contains eigenvectors associated with small eigenvalues.
In this section, we propose to build a subspace which does not excite the near-zero eigenvalues; by construction, the coercivity is preserved in that subspace and the condition number is controlled. We also prove that using this subspace only leads to a small loss of precision. By analogy with domain decomposition or multigrid methods, the subspace will often be referred to as the coarse space of approximation.

## 3.1 Notations

The analysis we develop is suited for discretizations in $L^2(C)$, either the Band discretization (15) or the Fourier discretization (16). An informal extension to the Dirac discretization (17) is proposed in section 3.5.

Note that the analysis can be conducted independently on each subdomain. For a given basis in the amplitude domain $\mathbf{A}^E$, we write $\boldsymbol{\alpha}^E$ the vector of components of the amplitude $A^E$, which is associated with the pressure field $p^E$:

$$ 
\begin{aligned}
s \in C &\mapsto A^E(s) = \mathbf{A}^E(s)\boldsymbol{\alpha}^E \\
x \in \Omega^E &\mapsto p^E(x) = \mathbf{P}_{\mathbf{A}}^E(x)\boldsymbol{\alpha}^E
\end{aligned}
$$

For the discretizations by piecewise constant amplitudes or by truncated Fourier series we note $\|\|_C$ the $L^2(C)$-norm of amplitude vectors and $\mathbf{M}_C^E$ mass matrix associated with basis $\mathbf{A}^E$:

$$ \|\boldsymbol{\alpha}^E\|_C^2 = \|\mathbf{A}^E\boldsymbol{\alpha}^E\|_{L^2(C)}^2 = \boldsymbol{\alpha}^{EH}\mathbf{M}_C^E\boldsymbol{\alpha}^E \quad (23) $$

Note that for any of the proposed discretization methods, the mass matrix $\mathbf{M}_C^E$ is diagonal.

## 3.2 Principle

We assume that for each subdomain $\Omega_E$, an approximation subspace was defined, and that, based on engineering rules presented in equation (18), it is sufficiently large for the solution to be approximated with enough accuracy.

We propose to select a smaller approximation subspace where amplitude generate non-negligible pressure. Note that this selection can be conducted independently on each subdomain, in parallel. For any subdomain $\Omega_E$, we would like to select the largest subspace $\tilde{\mathcal{A}}_\sigma^E$ of $\mathcal{A}^E$ where the following property is verified:

$$ \forall A \in \tilde{\mathcal{A}}_\sigma^E, \quad \|\mathcal{H}_E[A]\|_{S(\Omega_E)} \ge \sigma \|A\|_{L^2(C)} \quad (24) $$

where $\sigma > 0$ is a parameter to be defined by the user. In that subspace, the equation (22) becomes:

$$ \forall (A^E) \in \prod(\tilde{\mathcal{A}}_\sigma^E), \quad a_T((\mathcal{H}_E[A^E]), (\mathcal{H}_E[A^E])) \ge \sigma^2 Q \left( \sum_E \|A^E\|_{L^2(C)}^2 \right) \quad (25) $$

In other words, coercivity (and condition number) is controlled by $\sigma$. The smaller $\sigma$, the larger $\tilde{\mathcal{A}}_\sigma^E$ and the poorer the condition number is.

## 3.3 Construction of the subspace

The subspace $\tilde{\mathcal{A}}_\sigma^E$ is not practical to compute. We thus propose to a way to approximate it. The idea is to use the property of finite element interpolation to obtain a good estimation of the $S$ norm. Let us consider a mesh $T_h^E$ of $\Omega_E$. Let $(x_i^E)_{1 \le i \le N_h^E}$ denote the nodes of $T_h^E$ and $\psi^E = (..., \psi_i^E, ...)$ be associated matrix of shape functions (of degree $q$). $h$ stands for the maximal length of the edges of the mesh, it must be chosen in agreement with the characteristic length of the problem $\text{Re}(1/k)$. Anyhow this mesh is used to make interpolations and compute norms, not to approximate the solution so that it can be much coarser than recommended for finite element computations. Moreover the meshes $(T_h^E)_E$ do not need to be compatible at interfaces ($h$ could even be set independently on each subdomain).

Let $\mathbf{P}_{\mathbf{A},h}^E$ be the matrix of the pressure field generated by the basis $\mathbf{A}^E$ evaluated at the nodes of the mesh:

$$ \mathbf{P}_{\mathbf{A},h}^E = \left( \dots \quad \mathcal{H}_E[A_j^E](x_i^E) \quad \dots \right), \quad \text{row } i, \text{ column } j \quad (26) $$
A field of $S^{N_E}(\Omega_E)$ writes $\mathbf{P}_A^E \alpha^E$, its finite element interpolation is $\psi^E \mathbf{P}_{A,h}^E \alpha^E$, the classical properties of finite element interpolation ensure that the distance between these two fields can be controlled: there exist a constant $C^E$ which depends on the mesh and on $k$, and which can be decreased at will by decreasing the characteristic size $h$ of the mesh or increasing the degree $q$ of the interpolation [25], such that:

$$ \forall \alpha^E, \quad \|\mathbf{P}_A^E \alpha^E - \psi^E \mathbf{P}_{A,h}^E \alpha^E\|_{S(\Omega_E)} \le C^E \|\alpha^E\|_C \quad (27) $$

Let $M_h^E$ be the $S(\Omega_E)$-mass matrix associated with the mesh $T_h^E$, we have:

$$ \|\psi^E \mathbf{P}_{A,h}^E \alpha^E\|_{S(\Omega_E)}^2 = \alpha^{E^H} \mathbf{P}_{A,h}^{E^H} M_h^E \mathbf{P}_{A,h}^E \alpha^E \quad (28) $$

The approximation $\tilde{\mathcal{A}}_{\sigma,h}^E$ of $\tilde{\mathcal{A}}_\sigma^E$ defined in (24) is generated by all the combinations of amplitudes which create interpolated pressure fields of sufficient norm. In other words, we solve the following problem:

$$ \text{find } \alpha^E \in \mathbb{C}^{N_E} / \alpha^{E^H} \mathbf{P}_{A,h}^{E^H} M_h^E \mathbf{P}_{A,h}^E \alpha^E \ge \sigma^2 \alpha^{E^H} \mathbf{M}_C^E \alpha^E \quad (29) $$

This corresponds to the selection of eigenvectors associated with generalized eigenvalues larger than $\sigma^2$ for the system of hermitian positive definite matrices ($\mathbf{P}_{A,h}^{E^H} M_h^E \mathbf{P}_{A,h}^E, \mathbf{M}_C^E$).

Let $\tilde{\mathbf{V}}^E$ be the subset of $\mathbf{M}_C^E$-normalized eigenvectors $\mathbf{V}^E$ associated with eigenvalues larger than $\sigma^2$, $\mathbf{A}^E \tilde{\mathbf{V}}^E$ is a basis of $\tilde{\mathcal{A}}_{\sigma,h}^E$. If we assume that the mesh was built such that $\sigma > C^E$, we have:

$$ \begin{aligned} \|\mathbf{P}_A^E \tilde{\mathbf{V}}^E \tilde{\mathbf{\alpha}}^E\|_{S(\Omega_E)} &= \|\psi^E \mathbf{P}_{A,h}^E \tilde{\mathbf{V}}^E \tilde{\mathbf{\alpha}}^E + (\mathbf{P}_A^E \tilde{\mathbf{V}}^E \tilde{\mathbf{\alpha}}^E - \psi^E \mathbf{P}_{A,h}^E \tilde{\mathbf{V}}^E \tilde{\mathbf{\alpha}}^E)\|_{S(\Omega_E)} \\ &\ge (\sigma - C^E) \|\tilde{\mathbf{\alpha}}^E\|_C \end{aligned} \quad (30) $$

This inequality shows that the use of a mesh to estimate the $S$-norm does not prevent to control coercivity (25) as long as the mesh is fine enough for constant $C^E$ to be small with respect to $\sigma$. Indeed the coercivity relation becomes:

$$ \forall (\tilde{\mathbf{\alpha}}^E) \in \prod \mathbb{C}^{\tilde{N}_E}, \quad a_T((\mathbf{P}_A^E \tilde{\mathbf{V}}^E \tilde{\mathbf{\alpha}}^E), (\mathbf{P}_A^E \tilde{\mathbf{V}}^E \tilde{\mathbf{\alpha}}^E)) \ge (\sigma - \max_E C^E)^2 Q \left( \sum_E \|\tilde{\mathbf{\alpha}}^E\|_C^2 \right) \quad (31) $$

## 3.4 Loss of precision

Searching the solution in the subspace $\tilde{\mathcal{A}}_{\sigma,h}^E$ instead of $\mathcal{A}^E$, one expects a loss of precision. We thus need to estimate the ability of $\tilde{\mathcal{A}}_{\sigma,h}^E$ to approximate fields in $\mathcal{A}^E$. For any given components $\alpha^E$ of a vector in the approximation subspace $\mathcal{A}^E$, let $\tilde{\mathbf{\alpha}}^E = \tilde{\mathbf{V}}^E \tilde{\mathbf{V}}^{E^H} \alpha^E$ be its $\mathbf{M}_C^E$-orthogonal projection on the subspace $\tilde{\mathcal{A}}_{\sigma,h}^E$, we have:

$$ \begin{aligned} \|\mathbf{P}_A^E \alpha^E - \mathbf{P}_A^E \tilde{\mathbf{\alpha}}^E\|_{S(\Omega_E)} &\le \|\mathbf{P}_A^E \alpha^E - \psi^E \mathbf{P}_{A,h}^E \alpha^E\|_{S(\Omega_E)} \\ &\quad + \|\psi^E \mathbf{P}_{A,h}^E \alpha^E - \psi^E \mathbf{P}_{A,h}^E \tilde{\mathbf{\alpha}}^E\|_{S(\Omega_E)} \\ &\quad + \|(\psi^E \mathbf{P}_{A,h}^E \tilde{\mathbf{\alpha}}^E - \mathbf{P}_A^E \tilde{\mathbf{\alpha}}^E)\|_{S(\Omega_E)} \\ &\le C^E (\|\alpha^E\|_C + \|\tilde{\mathbf{\alpha}}^E\|_C) + \sigma \|\alpha^E\|_C \\ &\le (2C^E + \sigma) \|\alpha^E\|_C \le 3\sigma \|\alpha^E\|_C \end{aligned} \quad (32) $$

The extra error due to the truncation is thus directly controlled by $\sigma$ (for a mesh satisfying $\sigma > C^E$).

## 3.5 Informal extension for the discretization with Dirac's distributions

The discretization of the amplitude space by Dirac distributions (17) can not be conducted in $L^2(C)$. We extend the method to that case by arbitrarily setting $\mathbf{M}_C^E = \mathbf{I}$; the (improper) notation $\|A^E\|_{L^2(C)}$ then stands for the Euclidean norm of the components $\|\alpha^E\|_2$.
# 4 Practical considerations

## 4.1 Setting of the threshold $\sigma$

As seen earlier, the parameter $\sigma$ controls the coercivity (and then the condition number) and the attainable precision of the coarse problem set in $\tilde{\mathcal{A}}_{\sigma}^{E}$: the smaller $\sigma$ the larger the coarse problem, the poorer its condition number and the greater the precision of its solution.

We propose to choose $\sigma$ by an energy criterion: it is indeed possible to choose which fraction of the total energy is present in the coarse model.

For any of the choices of discretization, the basis $\mathbf{A}^E$, the mass-matrix for the amplitudes $\mathbf{M}_C^E$ is diagonal and easy to compute, we can at no extra cost transform the generalized eigenvalue system into a classical eigenvalue problem. The sum of the eigenvalues ($\theta_i$) is then the trace of the matrix.

The criterion is then defined by a scalar $0 < \beta \le 1$, so that $(1 - \beta)$ corresponds to the fraction of energy inside the coarse model:

$$ \sigma \text{ such that } \sum_{\theta_i > \sigma} (\theta_i) > (1 - \beta) \sum_{i=1}^{N_E} \theta_i = (1 - \beta) \text{trace}((\mathbf{M}_C^E)^{-\frac{1}{2}} \mathbf{P}_{\mathbf{A},h}^{E H} \mathbf{M}_h^E \mathbf{P}_{\mathbf{A},h}^{E} (\mathbf{M}_C^E)^{-\frac{1}{2}}) \quad (33) $$

In the examples, we used the following values for $\beta = 0.25, 0.1, 10^{-4}, 10^{-6}$.

## 4.2 Simplified evaluation of the $S$-norm

Regarding the evaluation of the $S(\Omega^E)$ norm, if the mesh $T_h^E$ is sufficiently fine and regular (no region is over or under meshed), then the condition number of $\mathbf{M}_h^E$ (for the Euclidean norm) is $O(1)$ [26], which means that the $\mathbf{M}_h^E$-norm is correctly approximated by the Euclidean norm. In that case, the generalized eigenvalue problem becomes the singular values decomposition of the matrix $\mathbf{P}_{\mathbf{A},h}^{E} (\mathbf{M}_C^E)^{-\frac{1}{2}}$. The sum of all singular values being equal to the (easy to compute) Frobenius norm of matrix $\mathbf{P}_{\mathbf{A},h}^{E} (\mathbf{M}_C^E)^{-\frac{1}{2}}$, a criterion similar to the one described in previous subsection can then be applied to define $\sigma$.

In the following examples, we use discretization in Dirac's (equation (17) and subsection 3.5) and the simplified evaluation of the $S$-norm, which is that case corresponds to computing the SVD of $\mathbf{P}_{\mathbf{A},h}^{E}$ (since $\mathbf{M}_C^E$ is set to identity). Unless stated otherwise, we use regular meshes with a characteristic size $h = \lambda/3$, where $\lambda$ is the wave-length.

## 4.3 Solver

The coarse subspace ($\tilde{\mathcal{A}}_{\sigma,h}^{E}$) can be directly used to find an approximate solution of the problem. Since the associated system is well conditioned and small, a direct solver can be employed, moreover we proved that the added error is limited. It can also be used as the coarse space of a multigrid method or of an augmented solver.

Let us briefly present how the subspace is used as an augmentation space for the LSQR Krylov solver [27] which is known to behave correctly on VTCR systems. We introduce the basis of the coarse subspace, obtained by the diagonal concatenation of the subdomains' basis:

$$ \tilde{\mathbf{V}} = \begin{pmatrix} \tilde{\mathbf{V}}^1 & 0 \\ 0 & \ddots \\ \mathbf{V}^{N_{\Omega}} \end{pmatrix} \quad (34) $$

LSQR solves for the normal equation:

$$ \mathbf{K}^H \mathbf{K} \alpha = \mathbf{K}^H \mathbf{f} \quad (35) $$

In that case, augmentation consists in ensuring that the residual $(\mathbf{f} - \mathbf{K}\alpha)$ is orthogonal to range$(\mathbf{K}\tilde{\mathbf{V}})$. This is classically implemented by an initialization/projection $(\alpha_0, \mathbf{X})$ method: one sets $\alpha = \mathbf{X}\hat{\alpha} + \alpha_0$ and solves the system in $\hat{\alpha}$ with a classical LSQR:

$$ (\mathbf{K}\mathbf{X})\hat{\alpha} = \mathbf{f} - \mathbf{K}\alpha_0 \quad (36) $$

where

$$ \alpha_0 = \tilde{\mathbf{V}}(\tilde{\mathbf{V}}^H \mathbf{K}^H \mathbf{K} \tilde{\mathbf{V}})^{-1} \tilde{\mathbf{V}}^H \mathbf{K}^H \mathbf{f} \\ \mathbf{X} = \mathbf{I} - \tilde{\mathbf{V}}(\tilde{\mathbf{V}}^H \mathbf{K}^H \mathbf{K} \tilde{\mathbf{V}})^{-1} \tilde{\mathbf{V}}^H \mathbf{K}^H \mathbf{K} \quad (37) $$
<figure><img src="image_3.png" /><figcaption>Figure 1: Scattering by a sound-hard cylinder (Section 5): (a) discretized computational domain, (b) real part of the exact solution at f=1500Hz, (c) real part of the exact solution at f=3500Hz.</figcaption></figure>

Note that the coarse matrix $(\tilde{\mathbf{V}}^H\mathbf{K}^H\mathbf{K}\tilde{\mathbf{V}})$ is coarse hermitian positive definite. Using this method, one expects a good initialization (which makes the initial residual small). Note that because the coarse problem controls the higher part of the spectrum whereas the lower part is very populated near zero, the rate of convergence can not be much improved by the augmentation.

In the following the direct solution in the coarse subspace is referred to as O-VTCR, the augmented LSQR is A-VTCR. Of course, O-VTCR corresponds to the initialization of A-VTCR.

# 5 Academic example: Scattering by a sound-hard cylinder

We first evaluate our method on an academic problem for which the analytical solution is known: the scattering of a plane wave by a sound-hard cylinder obstacle. We evaluate the classical VTCR equipped with various solvers: Matlab's direct solver "\" (which in that case corresponds to LU-solver), Moore-Penrose pseudo inverse (Matlab's "pinv"); and we compare them to the new reduced approach used either directly with Matlab's solver "\" (O-VTCR, which in that case corresponds to a Cholesky solver) or used as the augmented LSQR approach (A-VTCR), both for $\beta$ of equation (33) equal to $\beta = 0.25, 0.1, 10^{-4}, 10^{-6}$. We compare the characteristics of the systems (size, condition number, attainable precision) as well as the convergence of the iterative solvers.

The sound-hard cylinder obstacle $S_1$ has a radius $R_1 = 0.5$ m. The surrounding acoustic medium $\Omega$ is considered to be air ($\rho = 1.25 \text{ kg.m}^{-3}$, $c = 330 \text{ m.s}^{-1}$, and $\eta = 0$). It is truncated by a concentric circular surface $S_2$ of radius $R_2 = 6$ m. The problem (1) is solved with the following boundary conditions:

$$ \left| \begin{array}{l} \displaystyle \frac{\partial p}{\partial n} + \frac{\partial p_{scat}}{\partial n} = 0 \quad \text{over } S_1 \quad (a) \\ \\ \displaystyle p + i \frac{c}{\omega} \frac{\partial p}{\partial n} = 0 \quad \text{over } S_2 \quad (b) \end{array} \right. \qquad (38) $$

where $p_{scat}$ denotes a plane wave propagating in the direction $\theta = \pi$ and $\omega$ the circular frequency in rad.s$^{-1}$. Equation (38)(b) is an absorbing condition that approximates the Sommerfeld radiation condition.

The exact solution of the resulting problem is known to have the analytical form:

$$ p_{ex}(r, \theta) = \sum_{m=0}^{\infty} \frac{i^m}{2 - \delta_{0m}} \left( \frac{J_{m+1}(kR_1) - J_{m-1}(kR_1)}{H_{m+1}^{(2)}(kR_1) - H_{m-1}^{(2)}(kR_1)} \right) \cos(\theta) H_m^{(2)}(kr) \quad (39) $$

where $\delta_{0m}$ is the Kronecker's delta ($\delta_{0m} = 1$ if $m = 0$ and $\delta_{0m} = 0$ otherwise), $J_m$ is the order $m$ Bessel function of the first kind, and $H_m^{(2)}$ is the order $m$ Hankel function of second kind. Two different frequencies are considered f=1500 Hz and f=3500 Hz, corresponding respectively to 9 and 21 wavelengths in the diameter of the obstacle. The exact solutions are represented in Figure 1.

The domain $\Omega$ is decomposed into 18 sub-cavities as illustrated in Figure 1.
Errors $\epsilon[p]$ are evaluated using the following expression which incorporates contributions from pressure discontinuities across interfaces:

$$ \epsilon[p] := \frac{1}{\|p_{ex}\|_{S(\Omega)}} \left[ \|p - p_{ex}\|_{S(\Omega)}^2 + \sum_{\Gamma \in \mathcal{F}} \|[[p]]\|_{L^2(\Gamma)}^2 \right]^{1/2} \quad (40) $$

We use the discretization by uniform Dirac's distribution ("rays") (17). In order to study the influence of the initial refinement of the discretization of the amplitude space, the number of rays in the subdomains $N_E = \dim(\mathcal{A}^E)$, is set using the criterion (18) where the parameter $\mu$ varies from 0.1 (insufficient discretization) to 2 (a priori more than enough rays to represent the solution). The integrals needed in the calculation of the matrix $K$ are performed with a numerical quadrature using a point density corresponding to 30 points per wavelength.

## 5.1 Characteristics of the coarse model

Figure 2 (left) presents the evolution of the size of the problems and (right) the evolution of the condition number as functions of the underlying discretization $\mu$ of the amplitude space, at frequencies 1500 Hz (plain curves) and 3500 Hz (dashed curves).

The dimension of the approximation space $(\mathcal{A}^E)$ is linear in $\mu$ (see (18)). The optimized subspaces $(\tilde{\mathcal{A}}^E)$ almost coincide with $(\mathcal{A}^E)$ for small $\mu$ but they tend not to grow for $\mu > 0.75$: beyond this limit, increasing the size of the approximation space does not mean that more energy is present (in terms of pressure). This can in particular be observed by the small number of extra vectors required to capture a proportion of $(1 - 10^{-6})$ of the energy instead of $(1 - 10^{-4})$.

The condition number of the unfiltered system explodes around $\mu = 0.7$. This also corresponds to the existence of almost zero energy vectors in the approximation subspace. On the contrary, the optimized subspace experiences only a slight increase of the condition number which reaches at worse $10^4$.

<figure><img src="image_4.png" /><figcaption>Figure 2: Scattering by a sound-hard cylinder (Section 5): Evolution of the matrix size (left) and its condition number (right) with the parameter $\mu$ at the f=1500Hz (plain line) and f=3500Hz (dashed line) when using a direct VTCR approach (black) or the O-VTCR with $\beta = 0.25, 0.1, 10^{-4}$, and $10^{-6}$ (respectively orange, red, magenta and blue).</figcaption></figure>

## 5.2 Attainable precision by the coarse model

Figure 3 presents the residual, the pressure error (with respect to the analytical solution) and the norm of the solution, obtained by the solvers depending on the initial discretization (at frequencies 1500 Hz and 3500 Hz).
The dashed lines correspond to the use of a direct solver on the complete (black) or on the reduced models (O-VTCR for various $\beta$), the plain lines correspond either to a pseudo inverse on the full model or to an augmented LSQR-solver (A-VTCR). For that study the number of LSQR iterations is equal to the dimension of the search space; this choice is meant to illustrate the bad conditioning of the system since the attained residual is far from unit round-off. Setting the stopping criterion for augmented-LSQR is discussed in next subsection.

For insufficient discretization ($\mu < 0.7$), the coarse spaces practically coincide with the full space and condition number is low so that all solvers give comparable results (similar error in pressure, very small residual).

After a threshold (which depends on $\mu$ and $\beta$), the coarse spaces do not evolve, and the associated performance of reduced models (O-VTCR) does not improve: for instance, at 1500 Hz and for $\mu > 0.9$, $\beta = 10^{-4}$ gives a precision $\epsilon[p]$ of $5 \cdot 10^{-4}$ and $\beta = 10^{-6}$ gives a precision of $5 \cdot 10^{-5}$. Because of the better conditioning of the system, the precisions of the coarse models are much more stable than the precision of the full model with the same solver (black dashed curve).

We observe that the correlation between the residual and the error strongly depends on $\mu$; and that for sufficient discretizations, precision is often much better than the residual can tell (the residual is hardly different between $\beta = 10^{-1}$ and $\beta = 10^{-4}$ whereas the pressure error is 15 times smaller).

The norm of the amplitudes vector $\alpha$ explodes with direct solvers when the discretization is sufficiently fine for good precision to be attained, whereas the O-VTCR gives an amplitude vector of almost constant norm (independent of $\mu$).

If a reduced O-VTCR model gives insufficient precision, it can be improved either by choosing $\beta$ closer to zero or by doing augmented LSQR iterations (A-VTCR). The obtained precision is presented by plain curves. Whatever $\beta$, the final precision of the pressure field is comparable with the one obtained by the reference pseudo-inverse solver. Here again the residual is not an excellent indicator for the quality of the pressure solution. We observe that the solutions given by the coarse models have very small amplitudes and that doing iterations slightly increases that norm but anyhow the result remains much smaller than what is obtained by the Moore-Penrose pseudo-inverse which is supposed to give the minimal norm solution. In other words, the pseudo-inverse's slightly smaller residual is paid by a strong increase of the amplitudes.

## 5.3 Performance of the augmented solver

We now study the convergence of augmented-LSQR solvers (for various $\beta$) compared to classical LSQR. The underlying discretization is sufficiently fine for a solution of good quality be obtained ($\mu = 1.2$). Figure 4 shows the evolution of the norm of the residual and of the error in pressure during LSQR-iterations.

We observe that in terms of residual, curves are almost linear and parallel: augmentation does not increase the convergence rate but enables to start at a much smaller level of error. Note that in terms of pressure, the final precision of all methods practically coincides with the starting precision of A-LSQR with $\beta = 1e^{-6}$.

## 5.4 Sensitivity with respect to the mesh

In this subsection, we investigate the difficulties caused by the use of a mesh to interpolate the pressure and estimate its energy. The figure 5 how the properties of the coarse subspace evolves according to the mesh refinement. The characteristic length of the mesh $h$ is normalized by the wavelength $\lambda$. Note that for $\frac{h}{\lambda} > \frac{1}{2}$, it is known that the interpolation error is not controlled which means that we analyze the left part of the plot ($\frac{h}{\lambda} \le \frac{1}{2}$).

In the domain of interest, we observe (on the first row of Figure 5) that the dimension of the selected coarse space is almost constant (in particular for the finer $\beta$'s). More qualitatively, we measure the error in pressure from the resulting O-VTCR approach (second row of Figure 5). For $\beta \le 10^{-2}$, the error reaches a plateau when the mesh is fine enough $\frac{h}{\lambda} \le \frac{1}{3}$. The plateau is reached for finer meshes when $\beta = 10^{-1}$. Of course, the value of the plateau depends on $\beta$.

Finally we propose a measure of the evolution of the subspace itself. Let $\tilde{V}_r$ be a (Euclidean-orthonormal) basis of the reference coarse space computed with the finest mesh, and $\tilde{V}_h$ be a (Euclidean-orthonormal) basis of the coarse space computed with a mesh of dimension size $h$. Since dimension of the coarse space tends to increase with the refinement of the mesh, we propose the following measure for the distance between the spaces:

$$d = \|(I - \tilde{V}_h \tilde{V}_h^H) \tilde{V}_r\|_{Fro} \quad (41)$$

$d$ is thus the Frobenius norm of the projection of the reference space orthogonally to the current space. It somehow corresponds to the dimension of the supplementary subspace of $\text{span}(\tilde{V}_h)$ in $\text{span}(\tilde{V}_r)$. When normalized
by the size of the coarse subspace, we see that as long as $h \le \lambda$, the relative distance with the reference subspace is less than 1%.

This short study validates the choice $h = \frac{\lambda}{3}$ as a good instruction for the mesher, in particular when $\beta \le 10^{-2}$.

# 6 Numerical example: two dimensional car cavity

The VTCR, O-VTCR and A-VTCR are used to solve an acoustic problem for the car cavity depicted in Figure (Fig. 6). The cavity is filled with air ($\rho = 1.25 \text{ kg.m}^{-3}$, $c = 330 \text{ m.s}^{-1}$, and $\eta = 0$), it is excited by a uniform harmonic pressure in the front (boundary condition of type (1a) with $p_d = 1$). The front and rear windows are hard walls (boundary condition of type (1c) with $v_d = 0$), while an impedance condition of type (1b) is prescribed over all other boundaries (with $Z = i(\rho\omega)/1245$). The cavity is decomposed into 8 sub-cavities.

The problem is solved at three different frequencies 3000 Hz, 10000 Hz and 15000Hz using either the classic VTCR with a pseudo inverse solver, the O-VTCR and A-VTCR with $\beta = 10^{-1}$ and $10^{-6}$. The discretization is done using the approximation by Dirac's (17) associated with $\mu = 1.5$, the integrals needed to compute the matrix $K$ are performed numerically. The augmented LSQR has a stopping criterion of a relative residual of $10^{-3}$. The calculated pressure field and amplitude distribution of each sub-cavity are represented in the figures 7,8 and 9. The scales used for the amplitude distribution are given by the circles on the right of the drawing corresponding to the different methods.

<table><caption>Table 1: Two dimensional car cavity (Section 6): size and condition number of the algebraical system and relative error obtain with the different strategies</caption><thead><tr><th></th><th></th><th>5000 Hz</th><th>10000 Hz</th><th>15000 Hz</th></tr></thead><tbody><tr><td rowspan="2">VTCR</td><td>size</td><td>1020 × 1020</td><td>2020 × 2020</td><td>3022 × 3022</td></tr><tr><td>condition number</td><td>8.9 10<sup>18</sup></td><td>3.47 10<sup>19</sup></td><td>5.6 10<sup>19</sup></td></tr><tr><td rowspan="5">A/O-VTCR<br/>β = 10<sup>-1</sup></td><td>size</td><td>346 × 346</td><td>684 × 684</td><td>1022 × 1022</td></tr><tr><td>condition number</td><td>423</td><td>1.3 10<sup>3</sup></td><td>4.9 10<sup>3</sup></td></tr><tr><td>Relative error O-VTCR</td><td>0.305</td><td>0.2934</td><td>0.3311</td></tr><tr><td>Aug LSQR iterations</td><td>98</td><td>103</td><td>143</td></tr><tr><td>Relative error A-VTCR</td><td>0.0036</td><td>0.0629</td><td>0.0037</td></tr><tr><td rowspan="5">A/O-VTCR<br/>β = 10<sup>-6</sup></td><td>size</td><td>526 × 526</td><td>1008 × 1008</td><td>1483 × 1483</td></tr><tr><td>condition number</td><td>1.8 10<sup>3</sup></td><td>3.9 10<sup>3</sup></td><td>5.5 10<sup>3</sup></td></tr><tr><td>Relative error O-VTCR</td><td>0.0116</td><td>0.0479</td><td>0.0051</td></tr><tr><td>Aug LSQR iteration</td><td>9</td><td>6</td><td>18</td></tr><tr><td>Relative error A-VTCR</td><td>0.0012</td><td>0.0438</td><td>0.0032</td></tr></tbody></table>

One can see that the pressure fields are very similar between the VTCR and the A/O-VTCR, with an exception for the O-VTCR with $\beta = 0.1$. This visual impression is confirmed by the measure of the error given in the table 1 (the direct VTCR is considered to be the reference). However the distributions of amplitudes are different. The ones obtained with the pseudo inverse have extremely large norms and therefore the identification of the main directions of propagation is impossible. Both O-VTCR and A-VTCR lead to much smaller amplitudes and relatively clear main directions of propagation. The iterations of A-VTCR lead to noisier portraits of amplitudes than O-VTCR, especially at the highest frequency.

The table 1 also presents the size and condition number of the algebraical system obtained with the different approaches. The condition number of the original systems is extremely high which corresponds to almost zero eigenvalues whereas the condition number of the reduced system is controlled around $10^3$. As expected, for the O/A-VTCR, the precision is improved when reducing the value of $\beta$. One can observe that for the three considered frequencies, the A-VTCR leads to an accurate results even for large $\beta$.

The table 2 shows the computational time of the different steps of the resolution for the three methods at the three considered frequencies. For comparison, relative time with respect to classical direct VTCR is also provided. As one can see, the assembly of the VTCR matrix is relatively inexpensive in comparison to the resolution, making the resolution the bottleneck of the VTCR.

Note that our Matlab code is far from optimized. In particular the natural parallelism of the construction of the optimized space ($\hat{A}_{\sigma,h}^E$) is not exploited, and the augmented LSQR solver is crudely implemented. The CPU time for O/A-VTCR could thus easily be reduced.
<table><caption>Table 2: Two dimensional car cavity (Section 6): CPU time required for the different steps of the resolution for the VTCR, O-VTCR and A-VTCR.</caption><thead><tr><th></th><th></th><th>5000 Hz</th><th>10000 Hz</th><th>15000 Hz</th></tr></thead><tbody><tr><td rowspan="3">VTCR</td><td>Matrix Assembly</td><td>0.34 s</td><td>1.32 s</td><td>3.82 s</td></tr><tr><td>System solving</td><td>0.95 s</td><td>5.75 s</td><td>15.32 s</td></tr><tr><td><strong>Relative resolution time</strong></td><td><strong>1</strong></td><td><strong>1</strong></td><td><strong>1</strong></td></tr><tr><td rowspan="7">A/O-VTCR<br/>β = 10<sup>-1</sup></td><td>Construction Coarse space</td><td>0.14 s</td><td>1.05 s</td><td>2.02 s</td></tr><tr><td>Projection and Factorization</td><td>0.25 s</td><td>1.84 s</td><td>2.84 s</td></tr><tr><td>Total Resolution time O-VTCR</td><td>0.39 s</td><td>2.89 s</td><td>4.86 s</td></tr><tr><td><strong>Relative resolution time O-VTCR</strong></td><td><strong>0.41</strong></td><td><strong>0.50</strong></td><td><strong>0.31</strong></td></tr><tr><td>Aug LSQR time</td><td>0.98 s</td><td>3.57 s</td><td>5.49 s</td></tr><tr><td>Total Resolution time A-VTCR</td><td>1.37 s</td><td>6.46 s</td><td>10.37 s</td></tr><tr><td><strong>Relative resolution time A-VTCR</strong></td><td><strong>1.44</strong></td><td><strong>1.12</strong></td><td><strong>0.67</strong></td></tr><tr><td rowspan="7">A/O-VTCR<br/>β = 10<sup>-6</sup></td><td>Construction Coarse space</td><td>0.14 s</td><td>1.05 s</td><td>2.02 s</td></tr><tr><td>Projection and Factorization</td><td>0.43 s</td><td>2.74 s</td><td>5.54 s</td></tr><tr><td>Total Resolution time O-VTCR</td><td>0.57 s</td><td>3.79 s</td><td>7.76 s</td></tr><tr><td><strong>Relative resolution time O-VTCR</strong></td><td><strong>0.6</strong></td><td><strong>0.65</strong></td><td><strong>0.51</strong></td></tr><tr><td>Aug LSQR time</td><td>0.11 s</td><td>0.45 s</td><td>1.41 s</td></tr><tr><td>Total Resolution time A-VTCR</td><td>0.68 s</td><td>4.24 s</td><td>9.17 s</td></tr><tr><td><strong>Relative resolution time A-VTCR</strong></td><td><strong>0.71</strong></td><td><strong>0.73</strong></td><td><strong>0.59</strong></td></tr></tbody></table>

Yet we observe interesting performance for O-VTCR and A-VTCR. The speed-up of these approaches relative to the direct VTCR seems to increase for large problems. This is due to the fact that direct solvers have a cubic complexity with respect to the size of the problem whereas the iterative solvers' complexity is quadratic (when the number of iterations remains relatively small).

The sole cases where the new methods behave poorly, are when we try to achieve good precision with iterations starting with a poor coarse problem (A-VTCR with $\beta = 10^{-1}$) on small problems ($f \le 10000$ Hz).

In the end it seems that choosing a high precision coarse model is the most interesting strategy ($\beta = 10^{-6}$). In that case, the coarse problem is two times smaller than the original system, the solution in the optimized subspace (O-VTCR) is already of good quality (for a computational time of 49% of the reference time); if needed few iterations are sufficient to lower the residual (A-VTCR) at a limited cost (the gain in terms of time is then 41%).

# 7 Conclusion

When approximating Helmholtz equation, the Variational Theory of Complex Rays is a powerful alternative to the classical finite element because it is less subjected to the dispersion error. However the underlying representation of the unknown field in terms of amplitudes involves a compact operator which causes an accumulation of eigenvalues near zero, loss of coercivity and explosion of the condition number of the discrete system. In this paper we propose to filter the approximation space by creating a basis of amplitudes which have a significant contribution to the unknown in the appropriate norm. Note that the proposed method should adapt seamlessly to other methods which involve a compact operator in the representation of the unknown field.

In practice, a mesh is introduced in order to handle the pressure fields generated by distributions of amplitudes; this mesh is much coarser than would be required for a reliable finite element computation, typically the mesh used to visualize the solution is sufficient.

The selection process is associated with computing the highest part of the spectrum of a generalized eigenvalue problem independently on each subdomain. In most cases a simplification can be applied which leads to using singular value decomposition instead of generalized eigensystem. A truncation parameter must be introduced by the user, which controls the coercivity of the discrete sesquilinear form and the loss of precision with respect to the complete approximation space. This parameter is connected to the amount of energy represented in the filtered subspace. Numerical experiments prove that the criterion makes the selection process almost independent of the mesh (assuming it is reasonably refined with respect to the wavelength).

The filtered subspace can be used to obtain a good quality approximation of the solution with a system of reduced size and good conditioning. It can also be used as the coarse grid of an augmented Krylov solver. In
practice having 99.9999% of the energy inside the coarse space seems a good choice since both conditioning and
precision remain correct, if needed only few LSQR iterations will be required to achieve full precision.

Assessments proved that the method gave interesting performance in terms of precision, stability of the
amplitude solution and even CPU time (even with a non-optimized implementation).

**Acknowledgement:** The authors wish to thank Marc Bonnet and Martin Vorhalik for their helpful
discussions.

References

[1] O. C. Zienkiewicz, *The Finite Element Method*. McGraw-Hill, 1977.

[2] A. Deraemaeker, I. Babuska, and P. Bouillard, "Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions," *International Journal for Numerical Methods in Engineering*, vol. 46, pp. 471-499, 1999.

[3] A. Moiola and E. A. Spence, "Is the Helmholtz equation really sign-indefinite?," *SIAM Review*, vol. 56, no. 2, pp. 274-312, 2014.

[4] D. A. Di Pietro and A. Ern, *Mathematical Aspects of Discontinuous Galerkin Methods*, vol. 69 of *Mathématiques et Applications*. Springer, 2012.

[5] C. J. Gittelson and R. Hiptmair, "Dispersion analysis of plane wave discontinuous Galerkin methods," *International Journal for Numerical Methods in Engineering*, vol. 98, no. 5, pp. 313-323, 2014.

[6] O. Cessenat and B. Despres, "Application of an ultra weak variational formulation of elliptic PDEs to the two-dimensional Helmholtz problem," *SIAM Journal on Numerical Analysis*, vol. 35, pp. 255-299, 1998.

[7] C. Farhat, I. Harari, and L. Franca, "The discontinuous enrichment method," *Computer Methods in Applied Mechanics and Engineering*, vol. 190, pp. 6455-6479, 2001.

[8] W. Desmet, P. Sas, and D. Vandepitte, "An indirect Trefftz method for the steady-state dynamic analysis of coupled vibro-acoustic systems," *Computer Assisted Mechanics and Engineering Sciences*, vol. 8, pp. 271-288, 2001.

[9] P. Ladevèze, "A new computational approach for structure vibrations in the medium frequency range," *Comptes Rendus Académie des Sciences Paris, série II*, pp. 849-856, 1996.

[10] T. Strouboulis and R. Hidajat, "Partition of unity method for Helmholtz equation: q-convergence for plane-wave and wave-band local bases," *Applications of Mathematics*, vol. 51, pp. 181-204, 2006.

[11] H. Riou, P. Ladevèze, and B. Sourcis, "The multiscale VTCR approach applied to acoustics problems," *Journal of Computational Acoustics*, vol. 16, no. 4, pp. 487-505, 2008.

[12] L. Kovalevsky, P. Ladevèze, H. Riou, and M. Bonnet, "The variational theory of complex rays for three-dimensional Helmholtz problems," *Journal of Computational Acoustics*, vol. 20, 2012.

[13] T. Huttunen, J. Kaipio, and P. Monk, "An ultra-weak method for acoustic fluid-solid interaction," *Journal of Computational and Applied Mathematics*, vol. 213, no. 1, pp. 166-185, 2008.

[14] T. Huttunen, P. Monk, and J. P. Kaipio, "Computational aspects of the ultra-weak variational formulation," *Journal of Computational Physics*, vol. 182, no. 1, pp. 27-46, 2002.

[15] P. M. Morse, *Theoretical Acoustics*. Princeton university press, 1968.

[16] R. Dautray and J.-L. Lions, *Mathematical Analysis and Numerical Methods for Science and Technology*. Springer, 2000.

[17] J. T. Oden, I. Babuška, and C. E. Baumann, "A discontinuous hp finite element method for diffusion problems," *Journal of Computational Physics*, vol. 146, no. 2, pp. 491-519, 1998.

[18] D. Colton and R. Kress, "On the denseness of Herglotz wave functions and electromagnetic Herglotz pairs in Sobolev spaces," *Mathematical Methods in the Applied Sciences*, vol. 24, pp. 1289-1303, 2001.
[19] N. Weck, "Approximation by Herglotz wave functions," *Mathematical Methods in the Applied Sciences*, vol. 27, no. 2, pp. 155-162, 2004.

[20] L. Kovalevsky, P. Ladevèze, and H. Riou, "The Fourier version of the variational theory of complex rays for medium-frequency acoustics," *Computer Methods in Applied Mechanics and Engineering*, pp. 142-153, 2012.

[21] W. Desmet, B. van Hal, P. Sas, and D. Vandepitte, "A computationally efficient prediction technique for the steady-state dynamic analysis of coupled vibro-acoustic systems," *Advances in Engineering Software*, vol. 33, pp. 527-540, 2002.

[22] H. Riou, P. Ladevèze, B. Sourcis, B. Faverjon, and L. Kovalevsky, "An adaptive numerical strategy for the medium-frequency analysis of Helmholtz's problem," *Journal of Computational Acoustics*, vol. 20, 2012.

[23] B. Cockburn, J. Gopalakrishnan, and R. Lazarov, "Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems," *SIAM Journal on Numerical Analysis*, vol. 47, no. 2, pp. 1319-1365, 2009.

[24] N. Nguyen, J. Peraire, and B. Cockburn, "An implicit high-order hybridizable discontinuous Galerkin method for linear convection-diffusion equations," *Journal of Computational Physics*, vol. 228, no. 9, pp. 3232-3254, 2009.

[25] S. C. Brenner and R. Scott, *The Mathematical Theory of Finite Element Methods*, vol. 15 of *Texts in Applied Mathematics*. Springer, 2008.

[26] L. Kamenski, W. Huang, and H. Xu, "Conditioning of finite element equations with arbitrary anisotropic meshes," *Mathematics of Computation*, vol. 83, pp. 2187-2211, 2014.

[27] C. C. Paige and M. A. Saunders, "LSQR: An algorithm for sparse linear equations and sparse least squares," *ACM Transactions on Mathematical Software*, vol. 8, no. 1, pp. 43-71, 1982.
<figure><img src="image_5.png" /><figcaption>Figure 3: Scattering by a sound-hard cylinder (Section 5): Evolution of the residual, pressure error and solution's norm with the discretization refinement ($\mu$) at frequencies f=1500Hz (left) and f=3500Hz (right) for the direct VTCR approach (black) or the A-VTCR and O-VTCR with $\beta = 0.25, 0.1, 10^{-4}$ and $10^{-6}$ (respectively orange, red, magenta and blue).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Scattering by a sound-hard cylinder (Section 5): Evolution of the residue (left) and the error $\epsilon[p]$ (right) with the normalized number of iterations of the LSQR at the f=1500Hz (top) and f=3500Hz (bottom) for the direct VTCR approach (black) or the O-VTCR with $\beta = 0.25, 0.1, 10^{-4}$ and $10^{-6}$ (respectively orange, red, magenta and blue).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: Scattering by a sound-hard cylinder (Section 5): Dependence of the coarse space wrt the mesh (size, resulting precision, distance to reference coarse space) at frequencies 1500 Hz and 3500 Hz, and for various values of $\beta$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Geometry and boundary conditions of the two dimensional car cavity.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: Two dimensional car cavity (Section 6): pressure field and amplitude distribution at 5000Hz obtained with different resolution method: VTCR (top), O-VTCR (second line) and A-VTCR (third line), with $\beta = 10^{-1}$ (left) and $\beta = 10^{-6}$ (right)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: Two dimensional car cavity (Section 6): pressure field and amplitude distribution at 10000Hz obtained with different resolution method: VTCR (top), O-VTCR (second line) and A-VTCR (third line), with $\beta = 10^{-1}$ (left) and $\beta = 10^{-6}$ (right)</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9: Two dimensional car cavity (Section 6): pressure field and amplitude distribution at 15000Hz obtained with different resolution method: VTCR (top), O-VTCR (second line) and A-VTCR (third line), with $\beta = 10^{-1}$ (left) and $\beta = 10^{-6}$ (right)</figcaption></figure>