Rational-approximation-based model order reduction of Helmholtz frequency response problems with adaptive finite element snapshots

1.   Introduction

Many engineering applications, e.g., in structural dynamics, geophysics, seismology, acoustics or vibro-acoustics, require the numerical approximation of solutions to time-harmonic wave propagation problems over a range of frequencies. Due to oscillations in the analytical solutions, accurate numerical approximations of frequency responses are computationally expensive and time-consuming, already for moderate frequencies. Therefore, in the multi-query context, when multiple solves of the model have to be performed for different frequency values, the "naive" approach of discretizing the original problem as many times as needed is usually unaffordable.

Let us denote by $u(z)$ the exact solution to the time-harmonic problem under consideration for a given wavenumber $z \in \mathbb C$. Model order reduction (MOR) methods aim at alleviating the computational cost by producing an approximation of the solution map $z \mapsto u(z)$ or the frequency response map of a quantity of interest (QoI), usually given by some functional of the solution $z \mapsto y(z): = F(u(z))$. The produced approximation (the so-called surrogate) has to be close to the QoI and, at the same time, should be cheap to evaluate. Customarily, MOR methods rely on a two-phase procedure. The offline phase consists in the computation of a finite-dimensional basis of snapshots, entailing numerical solutions of the original boundary value problem for a set of frequency values (the sample points). Based on the offline information, the surrogate for the QoI is then assembled. These steps often require a considerable computational effort. However, they are performed only once, and the surrogate is then stored for later use during the online phase, where the surrogate is evaluated (in real-time) at any new frequency value of interest.

In the mathematical and engineering literature we can distinguish two families of MOR methods for frequency response problems. Projection-based techniques (see ^[6]), e.g., proper orthogonal decomposition (POD), the reduced basis method (RB), and the multi-moment matching method, are widely employed and extremely powerful; however, they require access to the original problem, particularly, to the stiffness and mass matrices as well as the forcing term. In contrast, the so-called non-intrusive MOR methods produce a surrogate relying only on the precomputed set of snapshots. We mention, e.g., the Loewner framework ^[30], Padé-based techniques ^[2,10], and minimal rational interpolation (MRI) ^[36,38]. Non-intrusive MOR methods usually display great flexibility, since, in principle, they allow one to construct a surrogate starting from snapshots obtained via a black-box solver (e.g., commercial closed-source software). The price to pay might be a reduced accuracy (compared to intrusive methods) for a fixed set of samples.

We underline that, for frequency response problems like the Helmholtz equation (which is the main focus of this work), the analytical solution $u(z)$ can be proven to be a meromorphic function of the complex wavenumber $z$. It is therefore sensible to look for its surrogate in the class of rational $\mathcal X$-valued maps (here $\mathcal X = H^1(\Omega)$ with $\Omega$ being the physical domain), as all the above-mentioned MOR approaches do. Moreover, if the QoI is a linear functional of the solution field, e.g., $y(z) = F(u(z))$ with $F \in \mathcal X'$, then it inherits a meromorphic structure from $u(z)$. The case of quadratic functionals (e.g., the energy-norm of the frequency response) is more involved, and needs to be treated separately.

In standard MOR techniques, the snapshots of the "truth model" are all computed on one discretization of the considered physical domain. In the specific framework we are handling, this might represent a big drawback. Indeed, the analytical solution of the Helmholtz equation oscillates more and more as the frequency increases, and it may exhibit local features, namely, a local resonance-type behavior, depending on the shape of the physical domain and the considered frequency values. Therefore, when a wide range of frequencies is considered, accuracy is guaranteed only by using a uniformly refined mesh, which must be sufficiently accurate for all frequency values of interest. This can potentially entail a waste of computational resources during the offline phase. Moreover, the a-priori identification of one grid yielding small FE errors for all samples of the wavenumbers is a challenging task.

In contrast, in spatially adaptive MOR approaches, each snapshot is taken on a mesh adapted to the local features at a given parameter (i.e., wavenumber) value, and belongs to a problem-adapted finite element (FE) space that contains the "optimal" number of degrees of freedom (DoFs) required to obtain a certain accuracy. Spatially adaptive MOR methods for coercive parametric partial differential equations (PDEs) have been the subject of recent mathematical research. In particular, we refer to ^[4], where the finite dimensional subspace for a greedy RB method consists of snapshots computed via an adaptive wavelet scheme. Adaptive FEs have been coupled with POD and greedy RB in ^[26,46] and ^[50], respectively.

The novelty of the present contribution resides in the introduction of spatially adaptive MOR methods for non-coercive elliptic boundary value problems, specifically, the parametric-in-frequency Helmholtz equation. The offline information, upon which the surrogate construction relies, is a set of snapshots computed by means of an $h$-adaptive FE method (here, $h > 0$ refers to the so-called mesh size, the maximum diameter of the elements of the FE mesh). The motivation that drives us to adaptive MOR methods is twofold. First, the use of $h$-adaptive FE snapshots, each living on a different mesh of the domain, allows one to save computational resources during the offline phase. Second, it presents an additional advantage: standard MOR methods aim at the approximation of the high-fidelity FE solution residing in an a-priori chosen mesh, not investigating how large the error between the high-fidelity solution and the analytical solution of the PDE is, but rather assuming such error to be small (uniformly over the parameter range). On the other hand, spatially adaptive MOR methods set the analytical solution as approximation target, bounding the FE error by means of a posteriori estimators. Notably, this framework allows one to balance the FE error (related to the choice of the mesh at each parameter value) and the MOR error (related to the number of snapshots and to the MOR strategy), resulting in an improved overall approximation efficiency.

On the other hand, the adaptive approach implies intrinsic difficulties: for instance, linear combinations of snapshots cannot be easily computed. In principle, to circumvent this issue, one could express all the snapshots as elements of some common FE space. However, for adaptive mesh refinement based on newest vertex bisection (NVB) ^[32,43], this would entail the construction of the so-called global mesh overlay ^[19,42], i.e., the smallest common refinement of all adapted meshes. In many cases, this entails a prohibitive computational effort and, more importantly, it goes against the main purpose of $h$-adaptivity. Therefore, in all the algorithms that we propose in this work, we strive to never construct the global mesh overlay. On the contrary, we will only require the evaluation of scalar products of pairs of snapshots, which is equivalent to building overlays of pairs of meshes. See also ^[26,46].

In principle, any numerical scheme well-suited for the Helmholtz equation at fixed accuracy (other than the $h$-adaptive FEM) can be employed offline during the sampling phase, and coupled with all the MOR techniques here presented. Typical examples are the $hp$-FEM, multiscale and extended schemes. However, from the algorithmical point of view, we are limited to numerical schemes for which the computation of scalar products of pairs of snapshots is feasible. This issue becomes particularly challenging when non-nested meshes or non-nested spaces are employed. Therefore, for simplicity, in the present work we choose to compute the snapshots by means of the piecewise affine adaptive FEM, for which we know how to compute the aforementioned scalar products.

When the QoI is a scalar, one possible way to avoid dealing with snapshots on different meshes is to construct a rational surrogate for the output itself by means of standard rational interpolation (SRI) methods. Off-the-shelf rational approximation methods like AAA ^[34], the Loewner framework ^[1,30], or vector fitting ^[27] can be applied to this aim. We summarize this class of approaches in Section 4.1.

On the other hand, the present contribution is dedicated to the development of MOR methodologies providing surrogates for function-valued QoI. With this aim:

● We extend the SRI approach to function-valued QoIs, particularly, traces and/or restrictions of the form $v(z) = u(z)|_\omega\in L^2(\omega)$ with $\omega$ being a part of the physical domain $\Omega$ or of its boundary $\partial\Omega$, or possibly even coinciding with the entire domain $\Omega$. We refer to the resulting method as $\mathcal V$-SRI. In the finite-dimensional "vector" case, this approach can be related to vector-valued rational approximation methods like set-valued AAA ^[33] or the MIMO Loewner framework ^[1]. Still, we focus here on the infinite-dimensional case, where the quantity to approximate is function-valued.

● We considerably improve the MRI method ^[36,38] by (ⅰ) formulating it in barycentric coordinates, which allow for enhanced numerical stability properties; (ⅱ) extending it to the $h$-adaptive FE setting. The resulting method is also referred to as $h$-adaptive MRI.

The rest of the paper is organized as follows. In Section 2, the parametric-in-frequency time-harmonic wave problems of interest for our discussion (interior and scattering problems) are introduced, and their meromorphic structure is highlighted. In Section 3, we recall the $h$-adaptive FE method (FEM). The core of the paper is Section 4. There, the $h$-adaptivity in the space variable is combined with the use of several rational-based MOR techniques. In Section 5, the $h$-adaptive POD is presented. In Section 6, we discuss the case in which the quantity of interest is a real-quadratic functional of the analytic solution. Finally, in Section 7, several numerical results are provided, to discuss the performance and highlight advantages and disadvantages of the considered approaches.

2.   Representative problems of interest

For a given complex wavenumber $k^2$, we look for $u(k^2)\in \mathcal X: = H^1_{\Gamma_D}(\Omega)$, the weak solution of the following interior Helmholtz boundary value problem

In Eq (2.1), $\Omega\subset \mathbb R^d$ ($d = 1, 2, 3$) denotes a bounded domain, whose boundary is partitioned into the (possibly empty) subsets $\Gamma_D$ and $\Gamma_N$. We denote by $\mathit{\boldsymbol{\nu}}$ the outward-pointing normal to $\Gamma_N$ and we assume that $f\in \mathcal X' = H^{-1}(\Omega)$ and $g_N\in H^{-1/2}(\Gamma_N)$. Problem (2.1) admits a unique weak solution $u(k^2)$ for all $k^2\in \mathbb C\setminus\Lambda$, $\Lambda$ being the set containing all (real, positive) eigenvalues of the Laplace operator with mixed Dirichlet-Neumann boundary conditions (see ^[10]).

Our ultimate target is the approximation of a quantity of interest (QoI) that depends on the weak solution $u(k^2)$, for all values of the (complex) parameter $k^2$ sweeping a given range of interest $Z$. In particular, let $F: \mathcal X\to \mathbb C$ represent a goal functional over $\mathcal X$. Then, the target quantity for fixed $k^2\in Z$ is the scalar

We will consider linear as well as real-quadratic functionals. For instance, we may take

or

(or weighted versions of them) with $\omega\subset\overline{\Omega}$ being, for example, an arbitrary subdomain or curve. Linear and real-quadratic functionals are quite often of interest in applications, representing, e.g., average displacements and vibrational energy, respectively.

For simplicity, we restrict our presentations to a $k$-independent forcing term $f$, $k$-independent Neumann datum $g_N$, and homogeneous Dirichlet boundary condition on $\Gamma_D$. If an inhomogeneous Dirichlet datum is given, by its lifting one falls back to a Helmholtz equation with homogeneous Dirichlet boundary conditions but with $k$-dependent right-hand side. As we showcase in our numerical examples, our discussion extends in a straightforward way to $k$-dependent forcing terms and Neumann data. However, one should note that, as the complexity of the forcing term and/or Neumann datum (particularly, with respect to $k$) increases, a surrogate of higher order – thus entailing a higher computational cost – is usually necessary to attain a prescribed approximation accuracy.

The interior Helmholtz problem Eq (2.1) can be extended to model exterior scattering problems. Indeed, let $D\subset \mathbb R^d$ ($d = 1, 2, 3$) denote a compact scatterer, whose surface has impedance $\zeta\in \mathbb C\cup\{\infty\}$ and inward-pointing normal $\mathit{\boldsymbol{\nu}}$, and define $\Omega = \mathbb R^d\setminus D$. For a given wavenumber $k\in \mathbb C$, the solution $u(k)\in \mathcal X = H^1(\Omega)$ of

models (in frequency domain) the wave scattered by $D$ when it is hit by the horizontal plane wave $e^{\iota kx_1}$. In practice, one can truncate $\Omega$: for instance, given $R > 0$ large enough, one can set $\Omega' = \Omega\cap B({\bf 0}, R)$, with $B({\bf 0}, R)$ denoting a ball of radius $R$ (see Figure 1), or $\Omega' = \Omega\cap[-R, R]^d$ (see Section 7.3). When writing the boundary value problem on the truncated domain $\Omega'$, the Sommerfeld radiation condition at infinity (the third equation in Eq (2.5)) is replaced by some approximation on the boundary $\partial\Omega'\setminus\partial D$, e.g., the first-order absorbing boundary condition $\partial_{\mathit{\boldsymbol{\nu}}}u(k) = \iota ku(k)$. The obtained variational problem admits a unique weak solution for all $k\in \mathbb C$, with ${\rm{Im}}(k)\geq 0$ (see ^[12]).

We note that, in the exterior case, we are writing $u(k)$, as opposed to the notation $u(k^2)$ that we used for the interior case, because the wavenumber $k$ appears also with the first power in Eq (2.5) (and also in the first-order absorbing boundary condition). This will allow us to use a common framework (with the argument $z$ of the solution map $u(z)$ denoting either $k^2$ or $k$) in describing the features of the solution map in the two cases (see Section 2.1).

To conclude this section, we underline that, while it is usually quite difficult to extend theoretical results from interior to exterior problems, our methods can be applied to scattering problems without the need for any modifications.

2.1.   Meromorphicity of the solution map

Given the parametric-in-frequency interior boundary value problem Eq (2.1), it is natural to introduce the solution map (or frequency response map) $u: \mathbb C\to \mathcal X$ that associates to each complex wavenumber squared $z = k^2$ the weak solution $u(z)$ of the corresponding Helmholtz boundary value problem. In ^[10] it was proven that the solution map is meromorphic over $\mathbb C$, i.e., $u$ is holomorphic over the whole complex plane, except at a set of isolated points $\Lambda = \{\lambda_i\}_{i = 1}^\infty$, where it displays pole-type singularities. In particular, $\Lambda$ is the set of eigenvalues of the (minus) Laplace operator on $\Omega$ with the considered boundary conditions. Moreover, each pole has order 1, regardless of its multiplicity as Laplace eigenvalue. More precisely, we have the following expansion:

with $\varphi_i\in \mathcal X$, $-\Delta\varphi_i = \lambda_i\varphi_i$, for all $i$, and $\{f_i\}_i\subset \mathbb C$ being coefficients depending on the forcing term and boundary conditions of the problem.

This result (except for the pole order) was extended to scattering problems Eq (2.5) in ^[12], where, in contrast, the solution map $u(z)$ is a function of the wavenumber, i.e., $z = k$. In this case, all the poles have strictly negative imaginary part. However, a characterization as precise as in formula (2.6) is not available: the poles might have order larger than 1, and the corresponding residues will not be $\mathcal X$-orthogonal. Still, the use of rational $\mathcal X$-valued surrogates is justified also in this setting.

Using a notation that encompasses both cases (interior and scattering problem), the frequency response map can be expressed as

where, for all $i$, $\lambda_i\in \mathbb C$ is a pole of $u$ with (finite) multiplicity $\mu_i$, and $\{\zeta_{i, \ell}\}_{\ell = 1}^{\mu_i}\subset \mathcal X$ are its corresponding generalized residues.

The meromorphicity of the frequency response map is inherited by QoIs $y$ that are linear functionals of $u$. The poles of $y$ are a subset of the poles of $u$. In particular, $\lambda_i$ is a pole of $y$ if at least one of its residues $\zeta_{i, \ell}$ satisfies $F(\zeta_{i, \ell})\neq 0$, i.e., if $\zeta_{i, \ell}$ is not in the kernel of $F$. Real-quadratic functionals, see, e.g., Eq (2.4), are more involved, and are discussed in Section 6. Due to the added difficulties described there, unless otherwise specified, from here onward we will assume that the QoI $y$ is a linear functional of the solution $u$.

Before proceeding, it is important to make the following observation.

Remark 2.1. Ritz- and Petrov-Galerkin projections of the problems above preserve the meromorphicity of the solution with respect to $z$. As a consequence, when, e.g., FEs are employed to discretize the PDE over a fixed mesh of $\Omega$, the FE discrete solution is meromorphic in $z$. However, the (countably) infinite analytic poles and eigenfunctions in Eq (2.6) are replaced by finite dimensional FE counterparts.

3.   Solution with $h$-adaptive FEM

In this section, we describe the $h$-adaptive FEM used to compute the snapshots that serve as inputs to the proposed MOR methods.

3.1.   Problem formulation and its FE discretization

Let us consider the boundary value problem

where $f \in L^2(\Omega)$, $g_N \in L^2(\Gamma_N)$ and $g_R \in L^2(\Gamma_R)$, with $\Gamma_D$, $\Gamma_N$, $\Gamma_R$ being relatively open and pairwise disjoint subsets of the boundary $\partial\Omega$ forming a partition of it, i.e., $\partial\Omega = \overline{\Gamma_D} \cup \overline{\Gamma_N} \cup \overline{\Gamma_R}$. This problem covers both Eq (2.1) (corresponding to the case $\Gamma_R = \emptyset$) and (2.5) with first-order absorbing boundary conditions (after truncation of the infinite exterior domain). The variational formulation of Eq (3.1) reads as follows: find $u \in \mathcal X = H^1_{\Gamma_D}(\Omega)$ such that

where the sesquilinear form $b: \mathcal X\times \mathcal X \to \mathbb C$ is given by

The variational problem is well-posed for $k\in \mathbb R$, if $\Gamma_R \neq \emptyset$ or if $\Gamma_R = \emptyset$ and $k^2$ is not an eigenvalue of the Laplace operator with mixed Dirichlet-Neumann boundary conditions.

To approximate $u$, we consider a standard $\mathcal X$-conforming FEM, namely, we perform the Galerkin projection on a finite dimensional subspace of $\mathcal X$ spanned by FE functions. (Note that this requires $\partial\Omega$ to be polygonal.) Let $\mathcal{T}_\bullet$ be a regular simplicial mesh of $\Omega$. We assume that the partition of $\partial\Omega$ is resolved by $\mathcal{T}_\bullet$ and denote by $\mathcal{E}_\bullet$ the set of facets (edges in 2D) of $\mathcal{T}_\bullet$. We consider the space of globally continuous and $\mathcal{T}_\bullet$-piecewise affine functions (P1-FEM), i.e.,

Then, the approximation $u_\bullet \in \mathcal X_\bullet$ of $u$ is the solution of the following finite-dimensional variational problem: find $u_\bullet \in \mathcal X_\bullet$ such that

It is well-known that Eq (3.3) can be ill-posed even if Eq (3.2) is well-posed. However, if Eq (3.2) is well-posed and $\mathcal{T}_\bullet$ is sufficiently fine, then Eq (3.3) admits a unique solution; see, e.g., [7, Proposition 1].

3.2.   A posteriori error estimation and adaptive algorithm

Error estimation techniques and mesh adaptive methods are powerful tools to accelerate the convergence of the FEM. Starting with the pioneering works ^[15,16], several error estimation strategies for the FEM have been proposed; see, e.g., the monograph ^[3]. For works focused on a posteriori error estimation for the Helmholtz equation, we refer, e.g., to ^{[8,9,17,21,44]}.

In this work, to estimate the error between $u$ and $u_\bullet$, we consider the classical residual-based a posteriori error estimator, i.e.,

where, for all $T \in \mathcal{T}_\bullet$, the computable local refinement indicators are given by

In Eq (3.4), $h_T = {\rm{diam}}(T)$, while $\left[\!\left[ {\partial_{\mathit{\boldsymbol{\nu}}} u_\bullet} \right]\!\right]$ denotes the jump of the normal derivative of $u_\bullet$ over the interface of two interior elements. If Eqs (3.2) and (3.3) are both well-posed, standard arguments in a posteriori error analysis (see, e.g., ^[3]HY__HY, Sections 2.2–2.3]) reveal that $\eta_\bullet$ is reliable and efficient in the sense that

Here, $C_{\mathrm{rel}}$ and $C_{\mathrm{eff}}$ are positive constants, while

collects the so-called data oscillations, where $f_T$ denotes the integral mean of $f$ in $T$, while $g_{N, e}$ (resp., $g_{R, e}$) denotes the integral mean of $g_N$ (resp., $g_R$) in $e \in \mathcal{E}_\bullet$. The positive constants $C_{\mathrm{rel}}$ and $C_{\mathrm{eff}}$ in Eq (3.5) depend on the shape-regularity of the mesh and on the problem data. They notably also depend on $k$.

Efficient error estimation is the fundamental ingredient to steer adaptive algorithms of the form

In such $h$-adaptive algorithms, the regions of the domain characterized by large values of the local error indicators are selected for further refinement. This mechanism automatically produces meshes that resolve the features of the solution and reduces the number of DoFs necessary to achieve a certain accuracy (compared to a uniform triangulation). For the Helmholtz equation, this results in meshes that resolve the local singularities of the analytical solution and are adjusted to the considered frequency value.

During the last twenty years, the (optimal) convergence of the adaptive FEM have been the subject of intense research; see, e.g., the recent review ^[18]. Such results guarantee that $h$-adaptive algorithms generate approximations with the desired accuracy in a finite number of iterations and that, asymptotically, the convergence has the best possible rate for the considered discretization. We refer to the recent work ^[7] for the proof of optimal convergence of an adaptive FEM for compactly perturbed elliptic problems such as the Helmholtz equation.

In the present work, for mesh refinement, we employ NVB ^[32,43]. Given a mesh $\mathcal{T}_\bullet$ and a subset $\mathcal{M}_\bullet \subseteq \mathcal{T}_\bullet$ of marked elements, we denote by $\mathcal{T}_\circ : = {\rm{refine}}(\mathcal{T}_\bullet, \mathcal{M}_\bullet)$ the coarsest NVB refinement of $\mathcal{T}_\bullet$ such that $\mathcal{M}_\bullet \subseteq \mathcal{T}_\bullet \setminus \mathcal{T}_\circ$. Furthermore, we assume that any mesh used to discretize $\Omega$ can be obtained by applying a finite number of NVB refinements to a given initial mesh $\mathcal{T}_0$.

Given the parameters $0 < \theta \le 1$, $\mathsf{tol}_h > 0$, and $N_{\mathrm{max}} > 1$, we consider the following algorithm (see [7, Algorithm 7]):

Apart from the outer if-statement, introduced in ^[7] to ensure that after a finite number of iterations the mesh becomes sufficiently fine to have well-posedness of Eq (3.3), Algorithm 1 is the standard adaptive algorithm of the form Eq (3.6). In step ($\star$), the selection of the elements to be marked for refinement, modulated by the parameter $0 < \theta \le 1$, is based on the Dörfler marking criterion Eq (3.7) from ^[20]. Algorithm 1 delivers an approximation $u_h = u_h(z) \approx u(z)$ such that either the corresponding error estimate $\eta_h$ satisfies $\eta_h \le \mathsf{tol}_h$ or the dimension of the underlying FE space satisfies $\dim \mathcal X_h > N_{\mathrm{max}}$.

3.3.   Computing bi-/sesquilinear forms across different meshes

The construction of the surrogate via the MOR methods discussed below will require the evaluation of scalar products of pairs of snapshots $u_h(z)$ and $u_h(z')$ associated with different frequencies $z$ and $z'$. This leads to some difficulties, since, in general, they belong to the different FE spaces $\mathcal X_h(z)$ and $\mathcal X_h(z')$. As such, it is necessary to assemble stiffness and mass matrices "across" different meshes $\mathcal{T}_h(z)$ and $\mathcal{T}_h(z')$. Exploiting the fact that $\mathcal{T}_h(z)$ and $\mathcal{T}_h(z')$ are NVB refinements of the same fixed initial mesh $\mathcal{T}_0$ and the binary tree structure of NVB, it is possible to compute the matrix entries efficiently and exactly, although the meshes are, in general, different and not even nested. For more details, we refer to [14, Section 6.2]; see also Section 7.

4.   Combining MOR with $h$-adaptivity

In the present section, we describe our proposed rational-based MOR methods for the approximation of $h$-adaptive quantities. First, we recall the SRI method in barycentric coordinates, which is suited for rational approximation of scalar-valued QoIs. Then, we extend the method to function-valued QoIs. Finally, we focus on the MRI method, and we specialize it to the adaptive framework of the present work.

4.1.   Standard rational interpolation (SRI) for scalar QoIs

When constructing the surrogate for a fixed QoI, a simple way to avoid dealing with snapshots on different meshes altogether is to approximate the (scalar) output $y_h(z) = F(u_h(z))$ directly. As observed in Section 2.1, linear functionals of $u$ are meromorphic, so that a rational approximation of $y_h$ is natural (see Section 4.4 for a further discussion on this):

with $P_{[N]}, Q_{[N]}\in\mathbb{P}_N(\mathbb C)$. Rational interpolation of scalar functions has been object of research for quite some time. A popular strategy entails the search of numerator and denominator of the rational approximant by minimization of the (weighted) linearized interpolation error at the sample points $\{z_j\}_{j = 1}^S\subset \mathbb C$, i.e.,

with suitable weights $\{w_j\}_{j = 1}^S\subset \mathbb R^+$. A normalization constraint (usually on $Q_{[N]}$) must be imposed to avoid the trivial solution $P_{[N]} = Q_{[N]} = 0$. A key property of rational interpolation (in this form) is that at least $2N+1$ samples are necessary. We note that, in the most general formulation of rational approximation, the numerator and denominator might have different polynomial degree. We ignore this here by allowing them to have defective degrees.

In the interest of numerical stability, a very useful representation of the polynomials $P_{[N]}$ and $Q_{[N]}$ can obtained in the so-called barycentric form:

where $\{\zeta_i\}_{i = 0}^N\subset \mathbb C$ are given distinct support points, and $\{p_i\}_{i = 0}^N\cup\{q_i\}_{i = 0}^N\subset \mathbb C$. Henceforth, we will refer to $\pi(z) = \prod_{i = 0}^N(z-\zeta_i)$ as nodal polynomial. Whenever a support point $\zeta_i$ coincides with a sample point $z_j$, interpolation at that point is guaranteed by setting $p_i = q_iy_h(\zeta_i)$, and the $j$-th addend can (and should) be removed from Eq (4.2). With this choice of basis, it is common to weigh the interpolation problem by nodal polynomial values: $w_j = |\pi(z_j)|^{-2}$. This allows to express Eq (4.2) in an extremely simple form.

Definition 4.1 (SRI). Let (distinct) sample points $\{z_j\}_{j = 1}^S$ and (distinct) support points $\{\zeta_i\}_{i = 0}^N$ be given, with $S\geq 2N+1$. The (barycentric) standard rational interpolant (SRI) of type $[N]$ of $y_h$ based on such sample and support points is a ratio $y_{h, [N]}^ {\rm{SRI}} = P_{[N]}^ {\rm{SRI}}\big/Q_{[N]}^ {\rm{SRI}}$, with $P_{[N]}^ {\rm{SRI}}$ and $Q_{[N]}^ {\rm{SRI}}$ of the form Eq (4.3). The coefficients $\{p_i\}_{i = 0}^N\cup\{q_i\}_{i = 0}^N\subset \mathbb C$ are chosen so that they minimize

under the constraints

(For more information on the barycentric rational form, see, e.g., ^[31,34], and the references therein.)

Similar (sometimes, equivalent) definitions yield some of the state-of-the-art algorithms for rational approximation, e.g., AAA ^[34], the Loewner framework ^[1], and vector fitting ^[27].

Remark 4.2. Let $i$ be such that $q_i\neq 0$. It is not difficult to see from Eq (4.3) that $y_{h, [N]}^ {\rm{SRI}}(\zeta_i) = p_i/q_i$. In particular, this means that, if $\zeta_i\in\{z_j\}_{j = 1}^S$ and $q_i\neq 0$, then we have interpolation of the target function at $\zeta_i$: $y_{h, [N]}^ {\rm{SRI}}(\zeta_i) = y_h(\zeta_i)$. This provides an intuitive justification for choosing the support points as a subset of the sample points $\{\zeta_i\}_{i = 0}^N\subset\{z_j\}_{j = 1}^S$, as done, e.g., in the AAA method ^[34].

Remark 4.3. In Definition 4.1, we have imposed a normalization on the Euclidean norm of the coefficients of the denominator $Q_{[N]}^ {\rm{SRI}}$ to exclude the trivial solution $P_{[N]}^ {\rm{SRI}} = Q_{[N]}^ {\rm{SRI}} = 0$. This choice has been observed to be more numerically robust than the standard one, based on the normalization of a single coefficient of the denominator, and is the de facto standard in rational approximation; see, e.g., ^[24,34].

Before proceeding, we outline here a practical algorithm for building SRIs. In the most common collocation case ($\{\zeta_i\}_{i = 0}^N\subset\{z_j\}_{j = 1}^S$), the coefficients $q_i$ can be found from the SVD of the Loewner matrix $G\in \mathbb C^{(S-N-1)\times(N+1)}$, defined entry-wise as

For more details, we refer either to ^[23,30,34], to Algorithm 2 in the next section, or to Remark B.2.

4.2.   SRI for function-valued quantities

The construction above can be easily generalized to vector quantities of interest, e.g., ${\bf y}_h(z) = \left(y^1_h(z), \ldots, y^O_h(z)\right)^\top$. To this aim, two changes are necessary: first, while the denominator remains $\mathbb C$-valued, the numerator must become $\mathbb C^O$-valued, i.e., $\{{\bf p}_i\}_{i = 0}^N\subset \mathbb C^O$; second, the absolute value in Eqs (4.2) and (4.4) must be replaced by the Euclidean norm over $\mathbb C^O$, yielding

We note that, in the collocation case ($\{\zeta_i\}_{i = 0}^N\subset\{z_j\}_{j = 1}^S$), this strategy is closely related (in fact, mostly equivalent) to the "set-valued AAA" and "fast-AAA" algorithms introduced in ^[33] and ^[29], respectively. An SVD-based solution is possible here as well, but we skip the details, since we provide them below in the more general framework of function-valued rational approximation, which includes this as special case. Alternatively, the MIMO Loewner framework ^[23,30] solves this kind of approximation problem by tangential interpolation, effectively recasting it as modified scalar problem.

More interesting and definitely less trivial is the extension to "infinite-dimensional" function-valued QoI. In the following, taking inspiration from Eqs (2.3) and (2.4), we focus on the case

However, as we describe in more detail in Section 4.4, we note that our discussion applies also to the case $\omega = \Omega$, i.e., $v = u$. We obtain the following definition by extending Eq (4.7) from $\mathbb C^O$ to $\mathcal{V}$.

Definition 4.4 ($\mathcal{V}$-SRI). Let (distinct) sample points $\{z_j\}_{j = 1}^S$ and (distinct) support points $\{\zeta_i\}_{i = 0}^N$ be given, with $S\geq 2N+1$. The (barycentric) $\mathcal{V}$-SRI of type $[N]$ of $v_h: \mathbb C\to\mathcal{V}$ based on such sample and support points is a ratio $v_{h, [N]}^{\mathcal{V} {\rm{-SRI}}} = P_{[N]}^{\mathcal{V} {\rm{-SRI}}}\big/Q_{[N]}^{\mathcal{V} {\rm{-SRI}}}$, with $P_{[N]}^{\mathcal{V} {\rm{-SRI}}}$ and $Q_{[N]}^{\mathcal{V} {\rm{-SRI}}}$ being of the form Eq (4.3). The coefficients $\{p_i\}_{i = 0}^N\subset\mathcal{V}$ and $\{q_i\}_{i = 0}^N\subset \mathbb C$ are chosen so that they minimize

under the constraints

Remark 4.5. In Eq (4.8), we have set $\mathcal V = L^2(\omega)$ with the objective of approximating some restriction of $u_h$ to $\omega\subset\Omega$. However, Definition 4.4 may be applied also in more general spaces $\mathcal{V}$. For instance, we may set $\mathcal{V} = \mathcal X$ (the space where the solution $u$ lives) to define $\mathcal X$-SRI.

It is crucial to observe that the following property (whose proof is deferred to Appendix A) holds true.

Lemma 4.6. The coefficients of the $\mathcal V$-SRI numerator $P_{[N]}^{ \mathcal V {\rm{-SRI}}}$ belong to the span of the snapshots: $\{p_i\}_{i = 0}^N\subset {\rm{span}}\{v_h(z_j)\}_{j = 1}^S$.

Accordingly, there exists a matrix $\mathring{P}\in \mathbb C^{S\times(N+1)}$ of expansion coefficients:

and the $\mathcal V$-SRI admits the expansion

An algorithm for computing $\mathcal V$-SRIs can be obtained as an extension of that for SRI. Notably, as we detail in Appendix B, the minimization of Eq (4.9) can be carried out through an SVD-like procedure in the "vectorized" $\mathcal V\otimes \mathbb C^S$ metric. As a practical way to do this, we introduce a $\mathcal V$-orthonormalization step: given $\{v_h(z_j)\}_{j = 1}^S\subset \mathcal V$, we compute a $\mathcal V$-orthonormal ($\langle\psi_j, \psi_{j'}\rangle_{ \mathcal V} = \delta_{jj'}$) basis $\{\psi_{j'}\}_{j' = 1}^T\subset \mathcal V$ and a $T\times S$ matrix $R$ such that

(The value $T\leq S$ is the rank of the snapshots, i.e., the dimension of their span as a subspace of $\mathcal V$.) The basis $\{\psi_{j'}\}_{j' = 1}^T$ can be constructed, e.g., by Householder triangularization of the snapshot quasi-matrix $[v_h(z_1)|\cdots|v_h(z_S)]$, see ^[45]. Alternatively, one can obtain $R$ from a Cholesky decomposition of a suitable Gramian matrix, cf. the proof of Lemma B.1.

We summarize the overall method in Algorithm 2. For simplicity, we only consider the framework where the support points are a subset of the sample points, with $\zeta_i = z_{i+1}$ for $i = 0, \ldots, N$. We note that the algorithm relies on the Loewner matrix $G\in \mathbb C^{T(S-N-1)\times(N+1)}$, defined entry-wise as

for $j = N+2, \ldots, S$, $j' = 1, \ldots, T$, and $i = 0, \ldots, N$, which is just a vectorized version of Eq (4.6).

The general case ($\{\zeta_i\}_{i = 0}^N\not\subset\{z_j\}_{j = 1}^S$) allows for a similar algorithm, where, however, the coefficients of the numerator are not necessarily scalar multiples of snapshots. We provide a strategy to compute $\mathcal V$-SRIs in Appendix B.

4.3.   Minimal rational interpolation (MRI) for state estimation

As we will motivate in Section 4.4, the approximation of the full state $u(z)$ is a topic that deserves particular attention per se. Due to its importance, several MOR techniques start by computing a surrogate for the state $u$ and only afterwards derive from it an approximation for the quantity of interest $y$. Among these, we can find projective approaches (see Section 5 and the references therein), least-squares Padé-based approaches (see ^[11]) and the previously mentioned $\mathcal X$-SRI. This last technique has been further developed in ^[38] under the name of minimal rational interpolation (MRI), with the objective of reducing the number of snapshots $S$ needed to achieve a certain rational type $[N]$ (or, vice versa, maximizing the rational type for a given number of snapshots). In particular, the MRI methodology allows constructing a rational approximant of type $[S-1]$, given only $S$ samples of $u$ (as opposed to the $2S-1$ that would be necessary in, e.g., SRI). In the following paragraphs, we not only extend MRI to the $h$-adaptive setting but we also formulate it in barycentric coordinates.

Definition 4.7 (MRI). The (barycentric) MRI of $u_h$ based on (distinct) sample points $\{z_j\}_{j = 1}^S$ is the ratio

$\pi$ being the nodal polynomial $\pi(z) = \prod_{j = 1}^S(z-z_j)$, where $\{q_i\}_{i = 0}^{S-1}\subset \mathbb C$ minimizes

under the constraint $\sum_{i = 0}^{S-1}\left|q_i\right|^2 = 1$.

Remark 4.8. As seen in the previous sections, the choice of coefficients of $P_{[S-1]}^ {\rm{MRI}}$ ensures interpolation of $u_h$ at all support points, which here coincide with all the sample points. Moreover, we observe that the target quantity Eq (4.15) corresponds to the $\mathcal X$-norm of the leading coefficient of the numerator:

It is interesting to note that the MRI approximant may be cast in the form Eq (4.12) used for $\mathcal X$-SRI (or $L^2(\Omega)$-SRI), where the coefficient matrix $\mathring{P}\in \mathbb C^{S\times S}$ is given by:

In fact, MRI could be interpreted as an extension of $\mathcal X$-SRI with an "unnatural" choice $S = N+1$.

The characterization of the denominator coefficients ${\bf q} = (q_0, \ldots, q_{S-1})^\top$ in the definition of MRI can be equivalently stated as: ${\bf q}$ is a minimal eigenvector of the snapshot Gramian

cf. Eq (B.1). This helps in designing a numerical strategy to compute the MRI surrogate, which we report in Algorithm 3.

4.4.   Some considerations on $\mathcal V$-SRI and MRI in the $h$-adaptive case

In the $h$-adaptive setting, samples of the quantities $v_h$ (see Section 4.2) at different parameter values live on potentially different meshes. The same holds for samples of the solution $u_h$. This represents an additional difficulty with respect to standard MOR (where the snapshots are all computed using the same mesh), which we must take it into account when defining, building, and evaluating the surrogate. For instance, evaluating the $L^2(\omega)$-SRI surrogate requires combining snapshots over all the (potentially different) meshes discretizing $\omega$, resulting in an overly high online cost due to a global mesh overlay. Unfortunately, this issue is intrinsic to the task of approximating a non-local quantity and cannot be easily mitigated. On the other hand, once $v_{h, [N]}^{L^2(\omega) {\rm{-SRI}}}\approx u|_\omega$ has been computed, it is easy to derive a posteriori (without training a new surrogate) surrogates for arbitrary linear functionals of $u(z)$, provided their support lies in $\omega$. One possible example is

with $w:\omega\to \mathbb C$ some weight function. In such cases, it suffices to extract the scalar samples $\{y_h(z_j)\}_{j = 1}^S$ from the infinite-dimensional ones $\{v_h(z_j)\}_{j = 1}^S$, and then replace^* $v_h$ by $y_h$ in Eq (4.12). Since the support of $F$ is contained in $\omega$, this can be done without using the full state $u_h(z)$.

^*Note that, in general, this yields a different rational surrogate from the one that would be obtained by SRI of $y$ directly. Indeed, the coefficients of the $L^2(\omega)$-SRI approximant (which we propose to reuse to obtain an approximation for $y$) are computed from information on $v$. On the other hand, the coefficients of the SRI surrogate of $y$ rely on information on $y$ only, which can be interpreted as a "filtered-down" version of $v$-information. This means that, while $L^2(\omega)$-SRI approximation followed by $y$-filtering is sure to be more versatile, approximating $y$ directly (by SRI) will likely be more accurate.

In view of this, computing the $\mathcal X$-SRI (or the MRI, or even the $L^2(\Omega)$-SRI) of $u(z)$ is quite powerful, since it allows to obtain a posteriori a surrogate for any linear functional of $u(z)$, without the need to repeat the training phase. As we will see in our numerical experiments in Section 7, the main price to pay for this flexibility is the offline computational cost needed to build the snapshot Gramian matrix Eq (B.1), which is necessary to find $R$. Indeed, each of its $\frac{N(N-1)}2$ independent non-diagonal entries requires building a (potentially different) pairwise mesh overlay over the whole $\Omega$, cf. Section 3.3.

Before proceeding further, we wish to mention that, while MRI has the advantage of using the samples optimally (and of being extremely simple to implement), it has one main drawback: its theoretical foundations hold only when the target of the surrogate admits a simple partial fraction decomposition whose residues are linearly independent. More explicitly, the following must hold true for the theory in ^[38] (particularly, the convergence of $u_{[S-1]}^ {\rm{MRI}}$ to $u$ in the $\mathcal X$-norm as $S\to\infty$) to apply:

$1)$ $u$ is of the form Eq (2.7), with simple poles, i.e., $\mu_i = 1$ for all $i$; note that the expansion has to hold only for $z$ in a (large enough) neighborhood of the target wavenumber range $Z$;

$2)$ the $S$ sample points $\{z_j\}_{j = 1}^S$ are judiciously placed on the wavenumber range of interest (good Lebesgue constant);

$3)$ the $N+1 = S$ most relevant residues (relevance is determined according to the Green's potential of $Z$, see ^[38]) should be linearly independent.

The first point is quite relevant to our discussion, since, in our setting, $u$ is of the required form, see Eq (2.6), but $u_h$ is not. A way to formalize this issue is the following: we take as approximation target the analytic solution $u$, whose snapshots are, however, affected by the noise $u_h-u$. From this viewpoint, an analysis of the properties of $h$-MRI could be based on the stability of MRI approximation with respect to noise. A theoretical discussion on this issue is planned to appear in a forthcoming paper. In the present work, we only observe the (sometimes adverse) effects of the $h$-FEM noise in our numerical experiments in Section 7. Of course, due to the adaptive nature of the snapshots, these noise-related issues are intrinsic to the problem at hand. As such, they may affect all the presented methods. Still, the LS-based formulations of SRI and $\mathcal V$-SRI have a regularizing effect that is missing in MRI.

5.   Projective MOR alternatives

In the literature, many MOR approaches have been proposed to approximate (functionals of) the solution of the Helmholtz equation. So-called projective methods build an approximation of the solution $u$, rather than of $y$ directly, through a restriction of the state problem (the Helmholtz equation) onto a suitably chosen subspace of $\mathcal X$, called "reduced space". Then the surrogate for $y$ is extracted from that of $u$, usually at small computational cost. The main feature characterizing a projective MOR technique is how such subspace is selected.

In POD, $u$ is evaluated at some prescribed parametric values, building a collection of snapshots. A principal components analysis of the snapshots is performed (e.g., by a generalized SVD), allowing the identification of their dominant modes. The reduced space is then defined as the span of such modes. We refer to ^[39] for an extensive discussion of POD. Then, the restriction of the original problem onto the reduced space is performed by (Petrov-)Galerkin projection onto $\widetilde{ \mathcal X} = {\rm{span}}\{\psi_j\}_j$, with $\{\psi_j\}_j$ the reduced basis of choice. Notably, this relies on a so-called affine decomposition of the original problem: Find $u(z)\in \mathcal X$ such that

gets projected onto $\widetilde{ \mathcal X}$ as: Find $\widetilde{u}(z) = \sum_ju_j(z)\psi_j\in\widetilde{ \mathcal X}$ such that

Obviously, this procedure requires access to each term $A_\ell$ and $b_\ell$ of the decomposition^†. Strategies are available to compute affine approximations of non-affine problems (see, e.g., ^[40]). In our $h$-adaptive setting, additional difficulties arise, since snapshots on different meshes are incompatible at the discrete level. Consequently, each of the operators $A_\ell$ appearing in Eq (5.1) must undergo a treatment similar to the one described in Section 3.3, see, e.g., ^[46].

^†As is evident from Eq (5.2), access to the (discretized) operator $A_\ell$ is usually necessary only through "matrix-vector multiplies", i.e., it must be possible to query $A_\ell\psi$ for an arbitrary $\psi\in \mathcal X$.

In the scope of this paper, we can highlight three main differences between MRI and POD:

● Approximation quality: for a fixed number of snapshots, Galerkin projection identifies a quasi-optimal approximation on their span, whereas for general problems no guarantee on the accuracy of MRI is available. In specific cases, e.g., for all Helmholtz problems of the form Eq (2.1), a theory for MRI can be derived, leading to quasi-optimality guarantees ^[38]. In such situations, the two approaches are often very comparable, cf. our numerical results in Section 7. In comparison, SRI or $\mathcal V$-SRI usually require twice as many snapshots to achieve a certain approximation order^‡.

^‡By "approximation order", for rational interpolation we mean the degree of the surrogate denominator plus one. For POD/RB, we mean the size of the reduced basis.

● Unfavorable problems: being based only on snapshots of $u$, MRI is somewhat oblivious to the structure of the underlying problem, which has $u$ as solution. As a consequence, some problems are inherently not amenable to approximation by MRI, while POD, being extremely structure-aware, manages to work well for them. This is an intrinsic difficulty for non-intrusive methods.

● Affinity requirements: the necessity for an affine (or affinely approximable by, e.g., EIM ^[39]) problem structure limits the applicability of POD. This is particularly relevant for scattering problems like Eq (2.5). Indeed, an affine decomposition of the impinging wave $e^{\iota kx_1}$ must necessarily include many terms, especially if the wavenumber range is large, thus hindering either online efficiency or surrogate accuracy (or both).

A schematic comparison of the presented algorithms, and of their computational costs, is included in Table 1.

6.   Rational interpolation for real-quadratic outputs

The case of real-quadratic functionals differs from that of linear ones because the meromorphic dependence on $z$ of the state is not inherited by real-quadratic outputs. Indeed, complex conjugation is not an analytic operation. For instance, consider the functional in Eq (2.4), applied to a solution map $u$ of the form Eq (2.7). We obtain

The quantity of interest $y(z)$ behaves like a meromorphic function for real wavenumbers only in special cases, namely, if the resonances $\lambda_j$ are real. Still, even in such situations, it is not meromorphic in the usual sense, since there is no complex neighborhood of the real axis where it is meromorphic.

In the more general case, we may ask whether it is possible to find a surrogate for the exact quadratic functional Eq (6.1) using any of the methods proposed above. The answer is actually positive for $L^2(\omega')$-SRI (provided $\omega\subset\omega'$), MRI, and (in many cases) POD. For the sake of simplicity, in the following paragraphs we restrict our attention to $L^2(\omega')$-SRI, with $\omega\subset\omega'$, but our discussion applies also to MRI.

Let Eq (6.1) be the quantity of interest, and let Eq (4.12) be the $L^2(\omega')$-SRI surrogate of $v_h = u_h|_{\omega'}$. Then, we observe that, by sesquilinearity,

where we have defined the $y$-representative Gramian $G_h^{[y]}$ as

In practical terms, this means that an online-efficient surrogate for $y$ can be built from the $L^2(\omega')$-SRI one, as long as the $y$-representative Gramian is assembled offline. Such assembly can be carried out by exploiting the strategy described in Section 3.3. Notably, we remark that $G_h^{[y]}$ can be built in conjunction to the $L^2(\omega')$-Gramian $G_h^{(v)}$ in Eq (B.1): in this case, the same overlay of two meshes can be used to compute $(G_h^{(v)})_{jj'}$, $(G_h^{(v)})_{j'j}$, $(G_h^{[y]})_{jj'}$, and $(G_h^{[y]})_{j'j}$, for $1\leq j < j'\leq S$.

Note that other similar approaches for handling real-quadratic outputs have been considered in "system theory" literature, e.g., (for projection-based MOR) in ^[47].

7.   Numerical results

In this section, we showcase the described methods in three numerical examples, where features and limitations of the different techniques emerge. We ran our code in MATLAB^® R2019b on a desktop computer with an 8-core 3.60 GHz Intel^® processor. The computation of the snapshots via the $h$-adaptive FEM is performed using $\texttt{p1afem}$ ^[22]. The assembly of stiffness and mass matrices associated with pairs of different meshes is carried out using [14, Section 6.2], with a computational complexity that is quadratic in the number of elements of the meshes. However, since NVB is a binary refinement rule, it should be possible to design an algorithm that performs this task in at most log-linear complexity with respect to the number of elements. This artificially inflates the computational times displayed below for the MRI and POD methods. Numerical testing with more efficient versions of the algorithm is planned for the upcoming future.

7.1.   Toy example on a triangle

We consider a triangular domain $\Omega = \left\{x\in \mathbb R^2, 0 < x_2 < x_1 < \frac{\pi}2\right\}$, whose boundary is partitioned into $\Gamma_1\cup\Gamma_2\cup\Gamma_3$. We are interested in the solution $u = u(z)\in H_{\Gamma_1}^1(\Omega) = \{v\in H^1(\Omega), v|_{\Gamma_1} = 0\}$ of the Helmholtz equation with uniform forcing term

The eigenproblem for the minus Laplacian over $\Omega$ with the above boundary conditions can be solved explicitly. The eigenvalues are

and the corresponding $L^2(\Omega)$-orthonormal eigenfunctions are

Accordingly, the solution $u$ admits the closed-form infinite series expression

where the series coefficients have been obtained by $L^2(\Omega)$-projection of the forcing term onto the eigenspaces. Hence, in the scope of evaluating the FEM and MOR errors in post-processing, the reference solution $u(z)$ is available to arbitrary precision by truncation of Eq (7.4). Below, we truncate the series by removing all terms in the sum whose coefficients $\frac{16}{\pi^2mn(m^2+n^2-z)}$ are smaller than $10^{-8}$.

7.1.1.   Adaptive FEM approximation for fixed frequency

In our first experiment, we present the computation of a snapshot using the $h$-adaptive FEM described in Section 3. To this end, we fix $z = k^2 = 51$ and run Algorithm 1 with $\theta = 0.1$, $\mathsf{tol}_h = 5\cdot 10^{-2}$, and $N_{\mathrm{max}} = \frac{\left|{\Omega}\right|k^4}{4 \mathsf{tol}_h^2}$. With this choice, we are heuristically trying to guarantee that the finest possible mesh used by the algorithm satisfies a "resolution condition" of the form $\widetilde{h}_\bullet k^2\sim 2 \mathsf{tol}_h$ ^[8], where $\widetilde{h}_\bullet^2 = \left|{\Omega}\right|/N_\bullet$ is the mesh size of a 2D uniform mesh made of right triangles. In this specific case, we have $N_{\mathrm{max}} \approx 3.2\cdot 10^6$, but the adaptive algorithm stops beforehand at $N_{143}\approx 1.9\cdot 10^5$, since $\mathsf{tol}_h$ is attained at the 143^rd iteration. Note that $z = 51$ is not of the form Eq (7.2), and its closest resonance is $\lambda_{1, 7} = \lambda_{5, 5} = 50$.

In Figure 2, we compare the analytical solution (Figure 2d) with the FE approximations at the 50^th (Figure 2a) and at the 143^rd iteration (Figure 2c). In Figure 2b, we show the evolution of the error estimator $\eta_\bullet$ and the true error $e_\bullet(z) = \left\|{\nabla(u_\bullet(z)-u(z))}\right\|_{{L^2(\Omega)}}$ as the mesh gets adaptively refined. Note that, in Algorithm 1, the $h$-adaptive loop stops as soon as $\eta_\bullet$ becomes smaller than the tolerance $\mathsf{tol}_h$. In Figure 2b, we are also including points that correspond to further adaptive steps of refinement only for illustrative reasons.

Several peaks appear before the asymptotic (optimal) convergence regime is reached. Such peaks are caused by resonances of the discrete problem at the chosen value of $z$. Indeed, since $z$ is not a resonance of the continuous problem Eq (7.1), we know that the discrete problem will not have a resonance at $z$ for a fine enough mesh. Still, resonances of the discrete problem may occur at $z$ if the mesh is too coarse. For accuracy, it is crucial that the adaptive algorithm stops after all peaks. We showcase this in Figure 2a, where we show the intermediate FEM solution after $50$ iterations, i.e., just to the left of the final peak in Figure 2b. We observe that the approximation of $u$ is rather poor, because the discrete resonance closest to $z$ is placed on the wrong side of $z$.

Instead, when the convergence enters the asymptotic regime, the error estimator and the true error decay with the same (optimal) rate. This is in agreement with the equivalence established in Eq (3.5) (note that here the oscillation terms vanish, because $f$ and $g_N$ are constant in this example).

Given $u$ and its FE approximation, we introduce the QoI

and the restriction $v(z): = u(z)\vert_{\Gamma_2}$ to $\Gamma_2$ (the domain of integration in Eq (7.5)). In Figure 3, we compare $v(z)$ with its FE approximation obtained by restricting $u^ {\rm{FEM}}_{51}(z)$ and $u^ {\rm{FEM}}_{143}(z)$ to $\Gamma_2$. We observe a complete mismatch between the exact solution and the coarse approximation, confirming our qualitative observations from the previous paragraphs. In Table 2, we make our conclusions quantitative by looking at the relative approximation error in the approximation of the QoI.

7.1.2.   Surrogate model for a linear quantity of interest

Now, we move to the approximation of the linear QoI $y(z)$ in Eq (7.5) over the range $z\in Z = [1, 100]$. First, we build a rational surrogate of type $[14/14]$ by SRI, using snapshots at $S = 29$ uniformly spaced sample points in $Z$. Specifically, we compute the snapshots by using the $h$-adaptive FEM with the same parameters ($\theta$, $\mathsf{tol}_h$, and $N_\text{max}$) as in the previous section. Note, in particular, that $N_\text{max}$ increases with $z$. When computing the snapshots at $z\approx 25.75$, we run into an issue: the $h$-adaptive loop in Algorithm 1 is terminated before the prescribed $\mathsf{tol}_h$ is attained, because the maximum number of elements is reached. This is caused by the nearby presence of the eigenvalue $\lambda_{1, 5} = 26$, which slows down the convergence of the FEM error. At this point, we are faced with a choice:

● We can use the inaccurate non-converged snapshot as it is. This is the simplest option, but might result in a poor surrogate, in particular near the affected sample point.

● We can discard the non-converged snapshot and build the surrogate without it. This is the "wasteful" option, since we are ignoring the results of valuable offline computations.

● We can (temporarily) increase $N_\text{max}$, e.g., multiplying its value by 10, and restart the $h$-adaptive iterations in the hope of convergence. If the increased $N_\text{max}$ leads to convergence, we use the snapshot. Otherwise, we discard it. This is the most robust, but potentially expensive, option. Here, we follow this approach.

Note that it is not necessarily a good idea to ditch $N_\text{max}$ altogether, or, equivalently, to set its value to $\infty$. Indeed, due to the a priori unknown location of the poles, it might happen that a sample point is extremely close to a pole of $u$. Correspondingly, the $h$-adaptive loop will necessarily take a long time to converge there. As such, we believe it more appropriate to be conservative, by setting a reasonably large (but finite) $N_\text{max}$, if necessary retrying with $10N_\text{max}$, and throwing away the snapshot in case of non-convergence.

In addition to the type $[14/14]$ SRI, we also build a type $[14/14]$ MRI and a type $[15/15]$ POD surrogate. Building the latter two reduced models requires only 15 snapshots, as opposed to the 29 needed for SRI. We take such snapshots at 15 uniformly spaced values of $z\in Z$. Note that, obviously, the issue of non-converging snapshots can appear also in MRI and POD. We deal with it as described above.

We show the results of the approximation in Figure 4. Note that the approximation error is computed with respect to the analytic solution Eq (7.4), so that, in particular, the error does not vanish at the sample points even though all approaches are interpolatory. This is due to the baseline FEM error, which acts as noise. We see that the approximations yielded by the three approaches are quite similar, with MRI behaving only slightly worse (on average) than the other two methods.

Looking locally around the (arbitrarily chosen) pole at $\lambda_{3, 5} = 34$, we see that SRI approximates best the location of the pole. This should be expected, due to the higher sampling resolution. Indeed, locally around $\lambda_{3, 5}$, SRI uses samples at both $z\approx 32.82$ and $z\approx 36.36$ while MRI and POD only have the one at $z\approx 36.36$.

In Table 3 we show the offline time of the three approaches. The computation of the snapshots in SRI takes considerably longer than for MRI/POD due to the non-converging samples discussed above. More generally, if the $S$ MRI/RB snapshots are a subset of the $(2S-1)$ SRI ones, the sampling for SRI will obviously take more time, but not necessarily twice as much. Overall, we see that building the SRI and MRI surrogates is about 20% faster than POD. Note that, if the SRI samples had all converged at the first try, we could have expected the offline phase of SRI to be, overall, about 4 times faster than that of MRI.

7.2.   Vibrations of a non-convex elastic plate

We consider a rectangular domain with a square hole

We denote by $\Gamma_1 = \left]0, 0.5\right[\times\{0\}$ and $\Gamma_2 = \left]0, 0.5\right[\times\{1\}$ the bottom and top sides of $\Omega$, respectively, see Figure 5. We are interested in the solution $u = u(z)\in H_{\Gamma_2}^1(\Omega) = \{v\in H^1(\Omega), v|_{\Gamma_2} = 0\}$ of the Helmholtz equation

The solution $u$ corresponds to the transverse displacement field of a thin membrane $\Omega$ in the "small deformation" regime, when the membrane, clamped at $\Gamma_2$, is subject to a time-harmonic excitation. The Neumann forcing term

denotes the force exerted on the membrane by the plane wave $u_\text{inc}(z)|_x = 0.1\exp\left(3\text{i}\sqrt{z}x\cdot\widehat{x}\right)$ (with $\widehat{x} = (\cos(\frac\pi6), \sin(\frac\pi6))^\top$ being the direction of propagation of the wave) that impinges on $\Gamma_1$ from below. Note that $g$ is not affine in $z$, cf. Eq (5.1), preventing a straightforward application of POD. We set the localized real-quadratic QoI

i.e., the mean-square displacement on $\Gamma_1$, as target of the approximation endeavor. This leads us to the definition of the intermediate quantity $v(z) = u(z)|_{\Gamma_1}\in L^2(\Gamma_1)$.

To obtain the FEM solution, we employ the $h$-adaptive FEM strategy described in Section 7.1.1. We show a sample FEM solution in Figure 5. For the approximation with respect to $z$, we employ $L^2(\Gamma_1)$-SRI and MRI, using $2S-1$ and $S$ uniformly spaced samples of $z\in[10, 200]$, respectively, with $S = 15$ and $S = 25$. This allows us to build rational approximations of the same type ($[S-1]$) with both methods.

The results are shown in Figure 6. There, we can observe that, for a given rational type, $L^2(\Gamma_1)$-SRI and MRI seem to achieve similar approximation errors. Notably, both approaches appear to struggle for low frequencies $z\approx 10$, especially in the cases where fewer snapshots are used.

A comparison of the corresponding runtimes can be found in Table 4. We can see that, for a given rational type, building the MRI method is more costly than the $L^2(\Gamma_1)$-SRI one, due to the computation of the snapshot Gramian. This is the case despite the lower amount of snapshots required by MRI. Note, in particular, that building the $L^2(\Gamma_1)$-Gramian is considerably faster than computing the $H_0^1(\Omega)$ one, due to the reduced dimensionality of the domain of integration.

7.3.   Acoustic scattering of a "trapping" cavity

We consider a domain of the form $\Omega = \Omega'\setminus D$, where $\Omega' = [0, 1]^2$ is a square and $D\subset\Omega'$ is a leftward-facing cavity, in the shape of a slanted "C", see Figure 7a. We are interested in the solution $u = u(z)\in H^1(\Omega)$ of the following Helmholtz equation with impedance boundary conditions (we set $z = k$ here):

The solution $u$ corresponds to the (time-harmonic) wave scattered by the sound-hard scatterer $D$ subject to a horizontal unit plane wave $u_\text{inc}(z) = e^{\iota zx_1}$ incoming from the left. The impedance condition on $\partial\Omega'$ serves as approximation of the Sommerfeld radiation condition at infinity, cf. Eq (2.5). Note that, as in the previous example, the Neumann datum is not affine in $z$. As approximation target, we take the trace $v(z) = u(z)|_\omega$, with $\omega\subset\partial D$ being the trapping boundary of $D$, i.e., the "interior" portion of the scatterer's boundary, see Figure 7b. Additionally, we consider the (squared) $L^2(\omega)$-norm of $v$, i.e., $y(z) = \int_\omega\left|{v(z)}\right|^2$.

For a given $z$, we compute $u(z)$ via the $h$-adaptive FEM strategy described in Section 7.1.1, but with $\mathsf{tol}_h = 0.5$. We show a sample FEM solution in Figure 7. As in the previous example, we employ $L^2(\omega)$-SRI and MRI for the approximation of $v(z)$ (and consequently $y(z)$) as $z$ varies in the interval of interest $Z = [10, 30]$. For both methods, we use $19$ uniformly spaced snapshots, allowing for the construction of rational approximations of types ^[9] and ^[18], respectively.

Before proceeding further, we note that, both with $L^2(\omega)$-SRI and with MRI, the approximation of $v(z)\in L^2(\omega)$ will require linear combinations of the snapshot traces $\{v(z_j)\}_{j = 1}^{19}$, which live on different discretizations of $\omega$. Still, since $\omega$ is a 1D curve, it is not too expensive to extend all the snapshot traces onto a common discretization of $\omega$, whose vertices are obtained as the union of the locations of all the DoFs of $\{v(z_j)\}_{j = 1}^{19}$. See also the discussion in Sections 4.2 and 4.3.

In Figure 8, we display the approximation of the scalar real-quadratic QoI $y(z)$ over the range of frequencies $z\in Z$. While the $L^2(\omega)$-SRI surrogate approximates the exact $y$ fairly well, we note the appearance of some oscillations in the MRI approximation, especially for $z\approx 30$. This behaviour is due to two compound effects: on one hand, MRI mistakenly places a surrogate pole quite close to the real axis, at $z\approx 29.1-0.14\iota$. On the other hand, modest Runge oscillations can be observed throughout the wavenumber range. We believe that the observed instabilities are due to the "interpolation" nature of the $h$-MRI method. Indeed, the employed offline information is a set of snapshots computed via the $h$-adaptive FEM, which are affected by the $h$-adaptivity-induced numerical noise $(u_h-u)$. Since the MRI surrogate interpolates the (noisy) snapshots, it is impacted by numerical instabilities. In contrast, the $L^2(\omega)$-SRI relies on a least-squares formulation, which helps in filtering out (part of) the numerical noise due to $h$-adaptivity.

In Figure 9, we show the approximation of $v(z')$ at the arbitrarily chosen frequency $z' = 25$, which is not one of the wavenumbers that are sampled for the construction of the surrogates. For the two rational surrogates, $\widetilde{v}$ is obtained simply by evaluating the rational approximations at $z'$. Looking at the error, we see that the $L^2(\omega)$-SRI surrogate is about one order of magnitude more accurate than the MRI one.

We show in Table 5 the timings corresponding to the two approaches. Our observations for the previous example hold true also here: overall, due to the necessity to compute the snapshot Gramian, building the approximant is more costly for MRI than for $L^2(\omega)$-SRI.

8.   Conclusions

In the present paper, we have presented several rational-based spatially adaptive MOR methods for the (non-coercive) parametric-in-frequency Helmholtz PDE endowed with mixed Dirichlet/Neumann/Robin boundary conditions. In particular, we have treated (ⅰ) SRI, which is useful in case of scalar-valued QoIs; (ⅱ) $\mathcal V$-SRI, which pertains to $\mathcal V$-valued QoI ($\mathcal V$ being a Hilbert space); (ⅲ) MRI, which provides a surrogate for the solution map $u(z)$ itself. The offline phase of each of the above-mentioned methods entails the solution of the considered problem for a set of values of the wavenumber. The snapshots are computed by means of the spatially adaptive FEM, and, as such, they reside in different discrete spaces, which are adapted to the value of the wavenumber and to the local features of each snapshot. As a projection-based alternative, we have also considered an $h$-adaptive version of POD.

With the target of comparing the methods, we have performed three numerical experiments. In the first example, we have observed that building the SRI and MRI surrogates is about $20\%$ faster than POD. The PDEs in the second and third example are both endowed with boundary conditions that depend non-affinely on $z$, preventing a straightforward application of POD. In the second example, the two methods are comparable, in the sense that they achieve similar approximation errors, with the $L^2(\omega)$-SRI method being faster than MRI. In our third example, the $L^2(\omega)$-SRI surrogate displays higher accuracy. Indeed, the MRI surrogate is affected by the presence of spurious poles and by numerical instabilities (Runge oscillations) due to the interplay between the two main features of the method, namely, $h$-adaptivity and interpolation of (noisy) snapshots. In contrast, the $L^2(\omega)$-SRI is more stable, since it produces an approximant in the least-square sense.

We remark that, in this work, the snapshots are computed via an $h$-adaptive algorithm designed to approximate the full state $u$ at optimal rate with respect to the energy norm. Another possibility, particularly appropriate if the functional representing the QoI is fixed, is to consider so-called goal-oriented adaptive algorithms, which aim at approximating at optimal rate the QoI directly; see, e.g., ^[25,35].

Moreover, all the algorithms are applied with sample points that are fixed a priori. In contrast, it is often desirable to follow a $z$-adaptive approach, where the selection of the sample points is driven by a suitable a posteriori estimator (note that this estimator would act on $z$, and not on the spatial domain, as $\eta_\bullet$ does). In the framework of projective MOR methods, this can be achieved with the weak-greedy RB technique. In some cases, the greedy selection of snapshots is also possible in the rational-based setting, e.g., in MRI, by means of the a posteriori estimator introduced in ^[36,37]. The derivation of novel MOR methodologies that combine $z$- and $h$-adaptivity is subject of a forthcoming work. In particular, with such an approach, we expect to be able to counteract some of the showcased instabilities (mainly, the spurious poles) of the $h$-MRI method.

An interesting further topic to be investigated is the extension of the presented $h$-adaptive MOR methods to the multi-parametric setting, i.e., Helmholtz-like problems where additional parameters (geometry, materials, etc.) are present on top of the wavenumber. This more general framework would be based on the parametric-MOR strategy presented in ^[36], and is of practical interest for several important applications, e.g., inverse problems, uncertainty quantification and optimal control problems, where non-$h$-adaptive multi-parametric MOR methods have been mostly employed so far (for instance, we mention ^{[5,13,28,48,49]}).

Acknowledgments

F. Bonizzoni is member of the INdAM Research group GNCS and acknowledges partial support from the European Research Council ERC under the European Union's Horizon 2020 research and innovation program (Grant agreement No. 865751). D. Pradovera acknowledges partial support from the Swiss National Science Foundation through project 182236. M. Ruggeri acknowledges partial support from the Austrian Science Fund (FWF) through the special research program Taming complexity in partial differential systems (grant F65). The authors would also like to acknowledge the kind hospitality of the Erwin Schrödinger International Institute for Mathematics and Physics, where part of this research was carried out.

Conflict of interest

The authors declare no conflict of interest.

A.   Proof of Lemma 4.6

In this section, we show that the coefficients of the $\mathcal V$-SRI numerator $P_{[N]}^{ \mathcal V {\rm{-SRI}}}$ belong to the span of the snapshots: $\{p_i\}_{i = 0}^N\subset {\rm{span}}\{v_h(z_j)\}_{j = 1}^S$.

Proof of Lemma 4.6. Without loss of generality, we assume that $N > 0$, and that sample and support points are sorted in such a way that there exists $s\in\{0, \ldots, N+1\}$ such that $z_j = \zeta_{j-1}$ for $j\leq s$, while $\{z_j\}_{j = s+1}^S\cap\{\zeta_i\}_{i = s}^N = \emptyset$.

For all $i = 0, \ldots, N$, consider a $\mathcal V$-orthogonal decomposition $p_i = p_i'+p_i"$, with $p_i'\in {\rm{span}}\{v_h(z_j)\}_{j = 1}^S$ and $\langle p_i", v_h(z_j)\rangle_{ \mathcal V} = 0$ for all $j = 1, \ldots, S$. First, we observe that, for all $i < s$, $p_i = q_iv_h(z_i)$ by Eq (4.10), i.e., $p_i" = 0$. So, if $s = N+1$, the claim is automatically proven. Otherwise, it remains to prove that $p_s" = \cdots = p_N" = 0$.

By the Pythagorean theorem, Eq (4.9) reads

For the sake of contradiction, assume that the optimal $p_s", \ldots, p_N"$ are not all equal to 0. By inspection of Eq (A.1), we see that setting all of them to zero cannot increase the value of the target, i.e., by optimality, the second sum must equal 0. However, this is absurd, since the Cauchy matrix

(which has more rows than columns) has full rank ^[41], so that

B.   General algorithm for the $\mathcal V$-SRI method

In Section 4.2 we have detailed the $\mathcal V$-SRI algorithm in the particular case where the support points are a subset of the sample points. A similar algorithm may be written even in the general case, where no relation between the set of sample points and support points is assumed.

Lemma B.1 ($\mathcal V$-SRI algorithm). Define the $v$-snapshot Gramian

whose rank is $T\leq S$. Also, let $s\in\{0, \ldots, N+1\}$ be the cardinality of $\{z_j\}_{j = 1}^S\cap\{\zeta_i\}_{i = 0}^N$. The corresponding $\mathcal V$-SRI exists and admits the following closed-form representation:

● The coefficients ${\bf q} = (q_0, \ldots, q_N)^\top$ of the $\mathcal V$-SRI denominator $Q_{[N]}^{ \mathcal V {\rm{-SRI}}}$ can be found as a (normalized) minimal right singular vector of a matrix $G\in \mathbb C^{T(S-s)\times(N+1)}$, i.e.,

● The coefficient matrix $\mathring{P}$, see Eq (4.11), may not be unique whenever $T < S$. On the other hand, one possible set of values for it can always be found as

In the proof, we provide expressions in closed form for $G$ and $H_i$, $i = 0, \ldots, N$.

Remark B.2. Lemma B.1 applies also to SRI, by replacing $v$ with $y$. Note that, in that case, $G_h^{(y)}$ has rank $T = 1$ and $G$ has size $(S-s)\times(N+1)$.

Proof of Lemma B.1. Without loss of generality, we assume that $N > 0$, and that sample and support points are sorted in such a way that $z_j = \zeta_{j-1}$ for $j\leq s$, while $\{z_j\}_{j = s+1}^S\cap\{\zeta_i\}_{i = s}^N = \emptyset$. Also, let $C$ and $\mathring{C}$ be the Cauchy matrices

and

respectively. By our assumptions on $N$ and $S$, $\mathring{C}$ has at least as many rows as columns, so that, in particular, it has full column rank ^[41].

For $i = 0, \ldots, s-1$, we have $p_i = q_iv_h(z_{i+1})$ by Eq. (4.10), so that, trivially,

and part of Eq (B.3) follows with

We now move to the cases $i = s, \ldots, N$.

Given $G_h^{(v)}$ in Eq (B.1), we first compute its rank-revealing Cholesky factorization $G_h^{(v)} = R^HR$, so that

Note that, given the matrix $R$, we can find a $\mathcal V$-orthonormal basis $\{\psi_{j'}\}_{j' = 1}^T\subset \mathcal V$ satisfying Eq (4.13). For instance, if $T = S$, we can set

(Note that such basis is never explicitly needed in the algorithm.)

In particular, by Eqs (4.11) and (4.13), we have the expansion

This, in turn, allows us to express Eq (4.9) in the Euclidean norm by the Pythagorean theorem

It suffices to take the gradient with respect to $\mathring{P}_{:i}$, for $i = s, \ldots, N$, to obtain the optimality conditions characterizing the numerator:

Using the previously introduced notation, these conditions can be equivalently expressed as

Denoting by ${\bf e}_{s+j}\in \mathbb C^S$ an element of the canonical basis ($({\bf e}_{s+j})_{j'} = \delta_{(s+j)j'}$), we have that $(G_h^{(v)})_{:(s+j)} = G_h^{(v)}{\bf e}_{s+j}$, and, for all $i = s, \ldots, N$,

This means that, component-wise,

for $i = s, \ldots, N$ and $j' = 1, \ldots, S$, with ${\bf z}_i = (z_{i1}, \ldots, z_{iS})^\top\in \mathbb C^S$ some arbitrary element of the kernel of $G_h^{(v)}$. We observe that the right-hand side above is affine in ${\bf q}$, so that we can write

with ${\bf a}_{ij'}\in \mathbb C^{N+1}$, defined entry-wise as

Now we collect the equations for a given $j'$ and $i = s, \ldots, N$:

By setting

we can write

This, with $\widetilde{{\bf z}}_i = {\bf 0}$, gives Eq (B.3) for $i = s, \ldots, N$.

Now, note that, since ${\bf z}_i$ belongs to the kernel of $G_h^{(v)} = R^HR$, ${\bf z}_i$ belongs to the kernel of $R$ as well. Thus, for all $i = s, \ldots, N$ and $j' = 1, \ldots, T$,

In particular, this means that, for a fixed ${\bf q}$, $\mathring{P}_{:i}$ may not be uniquely determined if $G_h^{(v)}$ is rank-deficient (i.e., if $T < S$), even though $p_i$ is always unique, see Eq (B.5).

Plugging Eq (B.8) into Eq (B.6), after proper re-indexing, yields

By setting

we deduce Eq (B.2), with