Localized Orthogonal Decomposition for two-scale Helmholtz-type problems

In this paper, we present a Localized Orthogonal Decomposition (LOD) in Petrov-Galerkin formulation for a two-scale Helmholtz-type problem. The two-scale problem is, for instance, motivated from the homogenization of the Helmholtz equation with high contrast, studied together with a corresponding multiscale method in (Ohlberger, Verf\"urth. A new Heterogeneous Multiscale Method for the Helmholtz equation with high contrast, arXiv:1605.03400, 2016). There, an unavoidable resolution condition on the mesh sizes in terms of the wave number has been observed, which is known as"pollution effect"in the finite element literature. Following ideas of (Gallistl, Peterseim. Comput. Methods Appl. Mech. Engrg. 295:1-17, 2015), we use standard finite element functions for the trial space, whereas the test functions are enriched by solutions of subscale problems (solved on a finer grid) on local patches. Provided that the oversampling parameter $m$, which indicates the size of the patches, is coupled logarithmically to the wave number, we obtain a quasi-optimal method under a reasonable resolution of a few degrees of freedom per wave length, thus overcoming the pollution effect. In the two-scale setting, the main challenges for the LOD lie in the coupling of the function spaces and in the periodic boundary conditions.


Introduction
The numerical solution of high frequency Helmholtz problems with standard finite element methods is still a very challenging task due to the oscillatory nature of the solutions. This manifests itself in the so-called pollution effect [3]: A much smaller mesh size than needed for a meaningful approximation of the solution are required for the stability and convergence of the numerical scheme. Typically, this leads to a resolution condition k α H = O(1) with α > 1 instead of kH = O(1), i.e. few degrees of freedom per wave length. The challenge becomes even greater when additionally studying the Helmholtz problem in a heterogeneous medium, such as locally periodic structures.
(Locally) periodic media, such as photonic crystals, can exhibit astonishing properties such as band gaps, artificial magnetism, or negative refraction [12,29,33,39]. One popular and successful modelling setup are scatterers made up of two materials with high contrast in the permittivity [4,5,6,7,10,16,28]. As the materials' fine-scale structures are much smaller than the wavelength, homogenization and corresponding numerical multiscale methods, such as the Heterogeneous Multiscale Method (HMM), are efficient tools to reduce the problem's complexity.
more sophisticated examples or even higher wave numbers may make the application of the LOD favorable, which hereby has the necessary theoretical footing. Concerning possible generalizations of our work, let us note that we concentrate on the two-scale setting described in [35], but the analysis can also be adapted to the case without high contrast [11] and might be useful for future numerical studies of the three-dimensional cases [4,6,7,10,28]. Moreover, the definition of the patches and interpolation operators in the periodic case and also the coupling and decoupling of various function spaces may be of general interest.
The paper is structured as follows: Section 2 introduces the notation and problem setting, see also [35]. In Section 3, the LOD in Petrov-Galerkin formulation is defined in the two-scale setting, where we also give some general comments on implementation aspects. Stability and convergence of the error are discussed in Section 4. Section 5 is devoted to the detailed proof of the decay of the correctors.

The two-scale Helmholtz problem
In this section, we introduce what we call the two-scale Helmholtz problem and the necessary notation. For further details on the derivation of this two-scale model and its practical relevance we refer to [35].
For the remainder of this article, let Ω ⊂⊂ G ⊂ R 2 be two domains with (polygonal) Lipschitz boundary at least. Throughout this paper, we use standard notation for Lebesgue and Sobolev spaces: For a domain ω, p ∈ [1, ∞), and s ∈ R ≥0 , L p (ω) denotes the complex Lebesgue space and H s (ω) denotes the complex (fractional) Sobolev space. The dot denotes a normal (real) scalar product, for a complex scalar product we explicitly conjugate the second component by using v * as the conjugate complex of v. The complex L 2 scalar product on a domain ω is abbreviated by (·, ·) ω and the domain is omitted if no confusion can arise. For v ∈ H 1 (ω), we frequently use the k-dependent norm v 1,k,ω := ( ∇v 2 ω + k 2 v 2 ω ) 1/2 , which is obviously equivalent to the H 1 norm. We write Y := [− 1 2 , 1 2 ) 2 to denote the twodimensional unit square. We indicate Y -periodic functions [2] with the subscript ♯. For example, can be interpreted as subspace of H 1 ♯ (Y ) and we write H 1 0 (D) ♯ to emphasize this periodic extension. By L p (Ω; X) we denote Bochner-Lebesgue spaces over the Banach space X.
The two-scale Helmholtz problem is now formulated as follows: We seek u := (u, with the two-scale function space H := H 1 (G) × L 2 (Ω; H 1 ♯,0 (Y * )) × L 2 (Ω; H 1 0 (D) ♯ ) and the twoscale sesquilinear form B defined by Beside the natural norm on H it is convenient for error estimation to define the following two-scale energy norm for a subdomain ω × R ⊂ G × Y with R 1 := R ∩ Y * and R 2 := R ∩ D. Furthermore, we introduce a version of the H 1 semi-norm on H via which is induced by the sesquilinear form We omit the subscript for the subdomain if ω × R = G × Y . In [35], we proved the following Lemma.
Lemma 2.1. There exist constants C B > 0 and C min := min{1, ε −1 e , Re(ε −1 i )} > 0 depending only on the parameters and the geometry, such that B is continuous with constant C B and fulfills a Gårding inequality with constant C min , i.e.
e . Note that C min can also be used to estimate the gradient terms in B by C min · 1,e from below. Furthermore, the unique solution u ∈ H to (2.1) (with additional volume term f ∈ L 2 (G) on the right-hand side) fulfills the following stability estimate see [35] for a proof with q = 3.
For the error analysis, we will compare the solution of the Localized Orthogonal Decomposition to a discrete reference solution, which is only needed for the theory and will never be computed in practical implementations. We introduce conforming and shape regular triangulations T H and T h of G and Y , respectively. Additionally, we assume that T H resolves the partition into Ω and G \ Ω and that T h resolves the partition into D and Y * and is periodic in the sense that it can be wrapped to a regular triangulation of the torus (without hanging nodes). We use the conforming We assume that this direct discretization is stable in the following sense: The (fine) mesh sizes H and h are small enough (in dependence on the wave number k) that there is a constant C HMM such that In [35], we discussed that this direct discretization can be re-cast into the traditional formulation of a Heterogeneous Multiscale Method and we proved that the discrete inf-sup-condition (2.6) holds under the classical resolution condition k q+2 (H + h) = O(1).

Remark 2.2.
We demonstrate the LOD at the specific example at the two-scale Helmholtz problem obtained in [35]. However, the theory can easily be extended to more general two-scale Helmholtz problems, which fulfill the following assumptions • the variational problem (2.1) involves a continuous sesquilinear form with Gårding inequality, i.e. an analogue of Lemma 2.1; • the analytical solution fulfills a stability estimate (2.4); • the (direct) Galerkin discretization (2.5) is stable (2.6).

The Localized Orthogonal Decomposition
In this section, we introduce the notation on meshes, finite element spaces, and (quasi)-interpolation operators and define the Localized Orthogonal decomposition in Petrov-Galerkin formulation for the two-scale setting. We close with some remarks regarding an implementation of the twoscale LOD.

Meshes and finite element spaces
Let the (fine) meshes T H of G and T h of Y be given as in the previous section, we assume that H and h are small enough that (2.6) is fulfilled. We consider a second, coarse discretization scale H c > H and h c > h: Let T Hc and T hc denote corresponding conforming, quasi-uniform, and shape regular triangulations of G and Y , respectively. As for the fine grids, we additionally assume that T hc is periodic and that T Hc and T hc resolve the partition of G into Ω and its complement and of Y into D and Y * , respectively. We denote by T hc (Y * ) and T hc (D) the parts of T hc triangulating Y * and D, respectively. The global mesh sizes are defined as H c := max{diam(T )|T ∈ T Hc } and h c := max{diam(S)|S ∈ T hc }. For the sake of simplicity we assume that T H and T h are derived from T Hc and T hc , respectively, by some regular, possibly non-uniform, mesh refinement including at least one global refinement. We consider simplicial partitions, but the theory of this paper carries over to quadrilateral partitions [17] and even meshless methods would be possible [22]. Given any subdomain ω ⊂ G define its neighborhood via N(ω) := int(∪{T ∈ T Hc |T ∩ ω = ∅}) and for any m ≥ 2 the patches see Figure 3.1 for an example. The shape regularity implies that there is a uniform bound C ol,m,G on the number of elements in the m-th order patch and the quasi-uniformity implies that C ol,m,G depends polynomially on m. We abbreviate C ol,G := C ol,1,G . The patches can also be defined in a similar way for a subdomain R ⊂ Y .
Here, we split R = R 1 ∪ R 2 with R 1 = R ∩ D and R 2 = R ∩ Y * , where R 1 or R 2 may be empty, and we write in short N m (R) := N m (R 1 )∪N m (R 2 ). N m (R 1 ) is defined in the same way as before, in particular, it ends at the boundary ∂D. For the patch N m (R 2 ) we interpret Y * as part of the torus. This implies that N m (R 2 ) ends at the inner boundary ∂D, but is continued periodically over the outer boundary ∂Y . This means that also the striped triangles in Figure 3.1 belong to the second patch for the periodic setting. We denote the overlap constants by C ol,m,Y and C ol,Y . By slight abuse of notation, we write We denote the conforming finite element triple space consisting of lowest order Lagrange elements with respect to the meshes T Hc and T hc by V Hc,hc as in the previous section. Again, we have V Hc,hc : ) and we moreover note that V Hc,hc ⊂ V H,h ⊂ H.

Quasi-interpolation
A key tool in the definition and the analysis is a bounded linear surjective (quasi)-interpolation operator I Hc,hc : V H,h → V Hc,hc that acts as a stable quasi-local projection in the following sense: It is a projection, i.e. I Hc,hc • I Hc,hc = I Hc,hc , and it is constructed as I Hc,hc := (I Hc , I Y * hc , I D hc ), where each (quasi)-interpolation operator fulfills the following. There exist con- Under the mesh condition that k(H c + h c ) 1, this implies stability in the two-scale energy norm The quasi-interpolation operator I Hc,hc is not unique: A different choice might lead to a different Localized Orthogonal Decomposition and this can even affect the practical performance of the method [37]. One popular choice is the concatenation of the l 2 projection onto piece-wise polynomials and the Oswald interpolation operator. Other choices are discussed in [14,38]. For the operators I Y * hc and I D hc not that they only act with respect to the second variable y. For I Y * hc , one can preserve periodicity as follows: The averaging process of the Oswald interpolation operator has to be continued over the periodic boundary (as for the patches before).

Definition of the LOD
The method approximates the discrete two-scale solution u H,h := (u H , u h,1 , u h,2 ) to (2.5) for given (fine) mesh sizes H, h. It is determined by the choice of the coarse mesh sizes H c and h c and the oversampling parameter m explained in the following. We assign to any (T, Note that the oversampling parameter does not have to be the same for G, Y * , and D. We could as well introduce patches N m0 (T ) × N m1 (S 1 ) × N m2 (S 2 ), but we choose m 0 = m 1 = m 2 =: m for simplicity of presentation and to improve readability.
We define the (truncated) finite element functions on the fine-scale meshes as where W H and W h are defined as the kernels of the corresponding (single) interpolation operators I HC and I Y * hc and I D hc , respectively. For given v Hc,hc ∈ V Hc,hc we define the localized correction The space of test functions then reads and can be written as triple We emphasize that dim V Hc,hc,m = dim V Hc,hc is low-dimensional and the dimension does not depend on H, h, or m.
The error analysis will show that the choice k(H c + h c ) 1 and m ≈ log k suffices to guarantee stability and quasi-optimality of the method, provided that the direct discretization (2.5) (with mesh widths H, h) is stable.
As discussed in [37], further stable variants of the method are possible: The local subscale correction procedure can be applied to only the test functions, only the ansatz functions, or both ansatz and test functions.

Remarks on implementation aspects
The present approach of the LOD exploits the two-scale structure of the underlying problem. In practice, one cannot work with the space triples such as V Hc,hc , but will look at each of the function spaces separately. The LOD consists of two main steps: First, the modified basis functions in V Hc,hc,m have to be determined, which includes the solution of the localized subscale corrector problems (3.3). Second, the actual LOD-approximation is computed as solution to (3.4). In this section, we explain how the computations in the macroscopic domain G and on the unit square Y can be decoupled in both steps. For general considerations on how to implement an LOD, for example algebraic realizations of the problems, we refer to [14].
Computation of modified bases. We observe that due to the sesquilinearity of B the following linearity for the correction operators Q m holds This means that the corrections of the basis functions in V 1 Hc , V 1 hc (Y * ) and V 1 hc (D) can be computed separately in the following way: Depending on the choice of the interpolation operator, Lagrange multipliers can be employed to decode that a function belongs to W H or W h , see [14].
We can further decrease the computational complexity of the localized corrector problems by decoupling the integrals over G and Y and by reducing the number of correction problems. The potential gain of course hinges on (additional) structure of the parameters and the meshes with the following general observations: • The corrections Q T ×S,m,1 and Q T ×S,m,2 only have to be computed for T ∈ T Hc with T ∩ Ω = ∅.
• In the case of constant parameters ε e and ε i , the corrector problems for Q T ×S,m,1 and Q T ×S,m,2 include information on T only in form of the weights |T | and |Ω T |; and the problems for Q T ×S,m only depend on S in form of the weights |S 1 | and |Y * S |. • In case of structured meshes T Hc and T hc and constant parameters, we can exploit symmetries to reduce the number of corrector problems [17].
Computation of the LOD-approximation. The LOD-approximation is defined as the solution to (3.4). This problem is similar to the discrete two-scale equation (2.5), only the test functions have been modified. Therefore, the LOD-approximation can be re-interpreted as an HMMapproximation with modified test functions and corrector problems. To be more explicit, u Hc,hc ∈ V Hc,hc from Definition 3.1 can be characterized as u Hc,hc = (u Hc , K hc,1 (u Hc ), K hc,2 (u Hc )), where u Hc ∈ V 1 Hc is the solution to a HMM with modified test functions and the corrections K hc,1 (u Hc ) and K hc,2 (u Hc ) are computed from u Hc and its reconstructions as described in [35]; see also [34,21] for similar reformulations in different settings. The HMM with modified test functions involves the following two steps: 1. Solve the cell problems for the reconstructions R 1 and R 2 around each quadrature point of the macroscopic triangulation T Hc using test functions in V hc,m (Y * ) and V hc,m (D).
2. Assemble the macroscopic sesquilinear form B H with the computed reconstructions and the test functions in V Hc,m .
Note that the reconstructions R 1 and R 2 as well as the fine-scale correctors K h,1 and K h,2 are different from those in [35] because of the modified test functions. This reformulation of the LOD-approximation as solution to a (modified) HMM decouples the computations on Y and G and no function triple spaces have to be considered. This is one great advantage of the present Petrov-Galerkin ansatz for the LOD in comparison to a Galerkin ansatz: We only need to compute reconstructions of standard Lagrange basis functions in V 1 Hc , but not of the basis functions in V Hc,m .

Error analysis
The error analysis is based on the observation that the localized subscale corrector problems The following result implies the well-posedness of the idealized corrector problems.  Proof. The essential observation is that for any (w h , w h,1 , w h,2 ) ∈ W H,h the property of the quasi-interpolation operators (3.1) implies that This directly yields the equivalence of norms on W H,h under the resolution condition (4.3). For the coercivity we observe that As the sesquilinear form B is also continuous (see Lemma 2.1), Lemma 4.1 implies that the idealized corrector problem (4.1) is well-posed and that the idealized correctors Q ∞ defined by (4.2) are continuous w.r.t. the two-scale energy norm  The stability of the LOD requires the coupling of the oversampling parameter m to the stability-/inf-sup-constant of the HMM. Therefore, we assume that H and h are small enough that (2.6) holds.  v Hc,hc e ψ Hc,hc e .
As C ol,m,G and C ol,m,Y grow at most polynomially with m because of the quasi-uniformity of T Hc and T hc , condition (4.6) is indeed satisfiable and the choice of the oversampling parameter m will be dominated by the logarithm of the wave number.
Proof. Let v Hc,hc ∈ V Hc,hc be given. From (2.6) we infer that there is ψ ∈ V H,h such that It follows from the structure of the sesquilinear form B that (Q ∞ (v * Hc,hc )) * solves the following adjoint corrector problem Hence, we obtain The second term on the right-hand side of (4.7) satisfies with Lemma 4.1 and Theorem 4.2 that |B(v Hc,hc , (Q ∞ − Q m )I Hc,hc ψ)| ≤ 1 + C min /2 C B (Q ∞ − Q m )I Hc,hc ψ 1,e v Hc,hc e ≤ 1 + C min /2 C B C 2 C ol,m,G + C ol,m,Y β m C I ψ e v Hc,hc e .
Hence, the condition (4.6) implies the assertion. The stability of the adjoint problem from Remark 4.4 implies z Hc,hc 1,e ≤ C LOD k q+1 C B I Hc,hc e e .
Thus, we obtain for (4.9) after division by I Hc,hc e e that I Hc,hc e e ≤ 1 + C min /2 C 2 B C 2 C ol,m,G + C ol,m,Y β m C LOD k q+1 · ( (1 − I Hc,hc )u H,h e + I Hc,hc e e ).
The oversampling condition (4.8) implies that the constants can be bounded by 1/2 and hence, the term I Hc,hc e e can be absorbed on the left-hand side. If the error for the direct (HMM-)approximation is small (which is the case for sufficiently small H, h), the error is dominated by the best approximation error of V Hc,hc , which can be quantified using standard interpolation operators and regularity results.

Proof of the decay property for the correctors
In this section, we give a proof of the exponential decay result of Theorem 4.2, which is central for this method. The idea of the proof is the same as in the previous proofs for the Helmholtz equation [8,17,37] or in the context of diffusion problems [22,30]. As in the previous sections, we have to take into account the two-scale nature of the problem and the spaces.
Let We note that periodicity is preserved when identifying degrees of freedom on the periodic boundary. Recall that the nodal Lagrange interpolation operator I is L 2 -and H 1 -stable on piece-wise polynomials on shape regular meshes due to inverse inequalities. Hence, for any (T, S 1 , S 2 ) ∈ T Hc × T hc (Y * ) × T hc (D) and all q ∈ P 2 (T ) × L 2 (T ; P 2 (S 1 )) × L 2 (T ; P 2 (S 2 )) we have the stability estimate I H,h q 1,e,T ×S ≤ C I q 1,e,T ×S . (5.1) In this section, we do not explicitly give the constants in the estimates. Instead we use a generic constant C, which is independent of the mesh sizes and the oversampling parameter, but may depend on the (quasi)-interpolation operators' norms, the overlap constants C ol,G and C ol,Y (not on C ol,m,G and C ol,m,Y !), the constant for the cut-off functions (see below), and C min .
In the proofs, we will frequently use cut-off functions. We collect some basic properties in the following lemma, cf. also [17, Appendix A, Lemma 2].
Since w := (1 − I Hc,hc )I H,h (ηφ) ∈ W H,h , the idealized corrector problem (4.1) and the fact that w has support only outside T × S imply B(w, φ) = B T ×S (w, v Hc,hc ) = 0. Therefore, we obtain Hence, estimates (5.5) and (5.4) give with the resolution condition (4.3) which gives the assertion after some algebraic manipulations.
As the localized correctors Q m are the Galerkin approximations of the idealized correctors Q ∞ , the decay property carries over to Q m . This is the main observation for the proof of Theorem 4.2.

Conclusion
In this paper, we presented a Localized Orthogonal Decomposition in Petrov-Galerkin formulation for the two-scale Helmholtz-type problems coming from homogenization, see for instance [35]. Under the natural assumption of a few degrees of freedom per wave length, k(H c + h c ) 1, and suitably chosen oversampling patches, m ≈ log(k), the method is stable and quasi-optimal without any further restrictions on the mesh width. Thereby, this work gives the theoretical foundation and justification for an application of the LOD to the two-scale setting in case the resolution condition poses a practical restriction for the direct discretization. We have demonstrated how the periodicity in the function spaces and the coupling of different spaces can be tackled in the LOD framework. Furthermore, this paper underlines the significance of viewing the Heterogeneous Multiscale Method as a direct discretization with numerical quadrature since only this viewpoint makes the additional application of the LOD possible.