Wavelet-based Edge Multiscale Finite Element Method for Helmholtz problems in perforated domains

TL;DR

This paper presents a wavelet-based multiscale finite element method for efficiently solving high-frequency Helmholtz problems in perforated domains, demonstrating convergence and effectiveness through numerical tests.

Contribution

It introduces a novel wavelet-based edge multiscale finite element method tailored for Helmholtz problems in perforated domains, allowing for large wavenumbers and providing proven convergence.

Findings

O(H) convergence under certain conditions

Effective for high wavenumbers in perforated domains

Validated by extensive 2D numerical experiments

Abstract

We introduce a new efficient algorithm for Helmholtz problems in perforated domains with the design of the scheme allowing for possibly large wavenumbers. Our method is based upon the Wavelet-based Edge Multiscale Finite Element Method (WEMsFEM) as proposed recently in [14]. For a regular coarse mesh with mesh size H, we establish O(H) convergence of this algorithm under the resolution assumption, and with the level parameter being sufficiently large. The performance of the algorithm is demonstrated by extensive 2-dimensional numerical tests including those motivated by photonic crystals.

Equations

$$\{\begin{array}{rlrl}-(\Delta +k^2)u& =f(x)& & \text{ in }{\Omega }^{ϵ}\\ \frac{\partial u}{\partial n}& =0& & \text{ on }\partial Q_1^{ϵ}\cap {\bar{\Omega }}_1\\ \frac{\partial u}{\partial n}-iku& =0& & \text{ on }\partial {\Omega }_2.\end{array}$$

$${\hat{Q}}_1^{ϵ}:=\{x:x/ϵ\in Q_1\}\text{and}{Q}_1^{ϵ}:={\hat{Q}}_1^{ϵ}+ϵ{\mathbb{Z}}^d.$$

$${\Omega }^{ϵ}:=({\Omega }_1\cap Q_1^{ϵ})\cup {\Omega }_2.$$

$${\Vert v\Vert }_V:=\sqrt{\Vert \nabla v{\Vert }_{{\Omega }^{ϵ}}^2+k^2\Vert v{\Vert }_{{\Omega }^{ϵ}}^2},$$

$$a(u,v)=(f,v{)}_{{\Omega }^{ϵ}}\text{for all }v\in V.$$

$$a(v_1,v_2):={\int }_{{\Omega }^{ϵ}}\nabla v_1\cdot \nabla {\bar{v}}_2\mathrm{d}x-k^2{\int }_{{\Omega }^{ϵ}}{v}_1\cdot {\bar{v}}_2\mathrm{d}x-ik{\int }_{\partial {\Omega }_2}{v}_1\cdot {\bar{v}}_2\mathrm{d}s\text{ for all }{v}_1,v_2\in V.$$

$$|a(v_1,v_2)|\le {\mathrm{C}}_{\mathrm{b}}\Vert v_1{\Vert }_V\Vert v_2{\Vert }_V\text{ for all }{v}_1,v_2\in V.$$

$$\mathrm{Re}\{a(v,v)\}+2k^2{\Vert v\Vert }_{L^2\left({\Omega }^{ϵ}\right)}^2\ge \Vert v{\Vert }_V^2\text{ for all }v\in V.$$

$${\Vert u\Vert }_V\le {\mathrm{C}}_{\mathrm{ap}}(k){\Vert f\Vert }_{L^2\left({\Omega }^{ϵ}\right)}.$$

$$a(v,z)=(v,w{)}_{{\Omega }^{ϵ}}\text{for all }v\in V,$$

$${\Vert z\Vert }_V\le {\mathrm{C}}_{\mathrm{ap}}(k){\Vert w\Vert }_{L^2\left({\Omega }^{ϵ}\right)}.$$

$$V_h:=\{v\in H^1({\Omega }^{ϵ}):V{|}_T\in {\mathcal{P}}_1(T)\text{ for all }T\in {\mathcal{T}}_h\},$$

$$a(u_h,v_h)=(f,v_h{)}_{{\Omega }^{ϵ}}\text{ for all }{v}_h\in V_h.$$

$${\omega }_i:=\bigcup \{K_j\in {\mathcal{T}}_H:\text{leavevmode}\text{leavevmode}\text{leavevmode}{O}_i\in {\overline{K}}_j\}.$$

$$C_{\mathrm{ov}}:=\underset{K\in {\mathcal{T}}_H}{\max }\#\{O_i:K\subset {\omega }_i\text{ for }i=1,2,\cdots ,N\}.$$

$$a(u_{\text{ms}},v)=(f,v{)}_{{\Omega }^{ϵ}}\text{for all}v\in V_{\text{ms}}.$$

$$x_{\ell ,j}=j\times h_{\ell },0\le j\le 2^{\ell }.$$

$$\varphi (x)=\{\begin{array}{rlrl}& 1,& & \text{ if }0\le x\le 1,\\ & 0,& & \text{ otherwise,}\end{array}\psi (x)=\{\begin{array}{rlrl}& 1,& & \text{ if }0\le x\le 1/2,\\ & -1,& & \text{ if }1/2<x\le 1,\\ & 0,& & \text{ otherwise}.\end{array}$$

$${\psi }_{\ell ,j}(x):=2^{\frac{\ell -1}{2}}\psi (2^{\ell -1}x-j)\text{ for all }0\le j\le 2^{\ell -1}-1.$$

$$W_{\ell }:=\{\begin{array}{rlrl}& \text{span}\{\varphi \}& & \text{ for }\ell =0\\ & \text{span}\{{\psi }_{\ell ,j}:0\le j\le 2^{\ell -1}-1\}& & \text{ for }\ell \ge 1.\end{array}$$

$$V_{\ell }:={\oplus }_{m\le \ell }{W}_m.$$

$$V_{\ell +1}=V_{\ell }{\oplus }_{L^2(I)}{W}_{\ell +1}$$

$$V_0\subset V_1\subset \cdots \subset V_{\ell }\subset V_{\ell +1}\cdots$$

$$L^2(I)={\oplus }_{\ell }{W}_{\ell }.$$

$$P_{\ell +1}v$$

$$\Vert v-P_{\ell }v{\Vert }_{L^2(I)}$$

$$-(\Delta +k^2){\chi }_i$$

$$\frac{\partial {\chi }_i}{\partial n}$$

$${\chi }_i$$

$${\mathcal{P}}_{i,\ell }(v):=\sum _{j=1}^{2^{\ell +2}}(v,{\psi }_j{)}_{\partial {\omega }_i}{\mathcal{L}}_i^{-1}({\psi }_j)\text{ for all }v\in L^2(\partial {\omega }_i).$$

$${\mathrm{C}}_{\mathrm{poin}}({\omega }_i)$$

$${\mathrm{C}}_{\mathrm{poin}}({\Omega }^{ϵ})$$

Full text

1 Introduction

The wave propagation through, and scattering from, complex multiscale structures is an important area of modern wave physics. The wave manipulation and control achievable by photonic crystals [18, 37], and more recent metamaterial devices [6, 8, 25, 35], underlie a host of wave devices in electromagnetism, optics and acoustics such as optical fibres, interferometers, mode converters, biosensors, thin-film optics for reflection control, optical switching and filtering and much more.

The canonical model problem is posed in terms of the Helmholtz equation in a perforated domain (see Figure 1):

[TABLE]

Here, we assume $f\in L^{2}(\Omega^{\epsilon})$, the wavenumber $k$ is real and positive, ${\mathbb{i}}$ is the imaginary unit, $\Omega^{\epsilon}$ is the perforated domain (the potentially complex structure as a model of a photonic crystal) to be defined in detail later, $\partial Q_{1}^{\epsilon}\cap\bar{\Omega}_{1}$ denotes the interface between the perforations and the perforated domain $\Omega^{\epsilon}$, and $\partial\Omega_{2}$ refers to the outer boundary. Complications arise in resolving the fine structure in the solution when $k$ is large, at high frequencies, and this is the regime often of interest in applications.

A key ingredient in any device design, or investigation of a physical effect, is the accurate, and fast, numerical simulation of the multiscale structures of interest. This has been the subject of intense and concentrated research over many years with a variety of techniques employed, the plane wave expansion methods [18, 19] are popular, as are multipole, Rayleigh, methods [37], finite difference time domain (FDTD) [20] and of course finite elements [37] figure strongly due to their versatility. Commercial finite element codes such as COMSOL [5] and FDTD such as Lumerical [24], dominate industry in terms of practicality, but there is a clear need for more modern implementations of finite-element based numerical methods in this field; many of the most interesting effects of topical interest are in three dimensions, such as flat lensing [9], antennas, or involve delicate changes in geometry, as in topological photonics [23], and may involve many tens, hundreds or even thousands of cells, each on the micro-scale, forming a macro-scale object where the wavelengths may be commensurate with the micro-scale; standard long-wave homogenisation is inappropriate and standard finite element approaches struggle to cope with the sheer size of memory required.

Separate from this area of physics, there has been extremely active research in the modern theory of finite elements and, in particular, on the development of efficient multiscale methods for practical applications with heterogeneous inseparable multiple scales. Due to this disparity of scales, classical numerical treatments become prohibitively expensive, and even intractable, for many multiscale applications. Nonetheless, motivated by the broad spectrum of practical applications, a large number of multiscale model reduction techniques, e.g., multiscale finite element methods (MsFEMs), heterogeneous multiscale methods (HMMs), variational multiscale methods, flux norm approach, generalized multiscale finite element methods (GMsFEMs) and localized orthogonal decomposition (LOD), have been proposed in the literature [2, 11, 10, 16, 17, 22, 27] over the last few decades. They have achieved great success in the efficient and accurate simulation of heterogeneous problems; we extend GMsFEMs using wavelets and furthermore investigate their application to the class of wave problems that encompass photonic crystals.

Designing efficient numerical solvers for Helmholtz equations with large wavenumbers has also attracted considerable attention over the past few decades; one of the main challenges is to reduce the so-called wavenumber dependent pollution effect [15, 30]. Mitigating the pollution effect, even for wave propagation through regular structures with homogeneous physical properties, requires an extremely fine mesh with grid size depending on the wavenumber $k$, or a very high polynomial degree $p$ within the basis. Consequently this results in an extremely expensive numerical scheme when the computational domain, or the wavenumber, is large. Numerical routes based even just around MsFEMs have not been explored in this context, and one would anticipate that they might result in efficient numerical solvers for wave propagation through complex multiscale structures; the WEMsFEM we develop fits broadly into the MsFEM framework but with generalisations and extensions. It is worthwhile noting that, for multiscale problems, the LOD approach has been investigated with [32, 33] having proposed numerical homogenization to eliminate the pollution effect for Helmholtz problems in heterogeneous media, however this is different from perforated domains and so their analysis does not carry over directly.

In this article we introduce a Wavelet-based Edge Multiscale Finite Element method (WEMsFEM) for Helmholtz equations in perforated domains, our Algorithm 1, inspired by the new multiscale algorithm proposed in [14, 21] for elliptic equations with heterogeneous coefficients. WEMsFEM takes advantage of the framework of GMsFEM [11], and utilizes the Partition of Unity Method (PUM) [36] as the essential component, and extends these approaches with additional novel ingredients that include a rather cheap local solver and provable convergence rate [14].

The main challenges in Problem (1.1) lie in accurately describing the interfaces, possibly a large number of perforations, and a large computational domain. The main idea of WEMsFEM is to utilize wavelets to approximate the solution restricted on the coarse edges, and then transfer this approximation property to the interior error estimate. Note that the coarse mesh cannot resolve the interfaces, nonetheless under the Scale Resolution Assumption (3.1) that the coarse mesh grid $H=\mathcal{O}(k^{-1})$, with $k$ being the wavenumber, we will prove in Proposition 4.1 that the error of our proposed multiscale algorithm in the energy norm is of $\mathcal{O}(H)$ given the wavelet parameter $\ell=\mathcal{O}(\log_{2}(kC_{\text{ap}}(k)))$ with $C_{\text{ap}}(k)$ being a stability constant defined as in (2.5).

The remainder of this paper is constructed as follows: We first present, in Section 2, the detailed problem and its basic properties. The Wavelet-based Edge Multiscale Finite Element Method (WEMsFEM) is introduced in Section 3 to solve this problem, and its convergence is analyzed in Section 4. Furthermore, we present in Section 5 extensive numerical tests to demonstrate the performance of our proposed algorithm. Finally, we draw together our results for discussion in Section 6.

2 General setting

In this section we present the general setting of the Helmholtz problem in a perforated domain; we also provide basic properties and results pertinent to the problem, and an outline of the construction of an ansatz space based on GMsFEM.

We start with the geometric setting of the domain for Problem (1.1). Let $Q:=[0,1]^{d}$ be the reference periodicity cell in $\mathbb{R}^{d}$ with $d\geq 2$ and we take $Q_{0}\subset Q$ with infinite smooth boundary $\partial Q_{0}$. Denoting $Q_{1}:=Q\backslash\bar{Q}_{0}$ as one unit cell with size $1$, then the contracted set $\widehat{Q}_{1}^{\epsilon}$ is one cell of the crystal with size $\epsilon$; $\widehat{Q}_{1}^{\epsilon}$, and its $\epsilon$-periodic cloning, $Q_{1}^{\epsilon}$, are defined as

[TABLE]

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded Lipschitz domain, and also let $\Omega_{1}\subset\Omega$ and denote by $\Omega_{2}:=\Omega\backslash\Omega_{1}$, then the computational domain is

[TABLE]

We introduce the complex-valued space $V:=H^{1}(\Omega^{\epsilon};\mathbb{C}):={W^{1,2}(\Omega^{\epsilon};\mathbb{C})}$, equipped with the $k$-weighted norm

[TABLE]

and similarly define the complex-valued space $V(D):=H^{1}(D;\mathbb{C})$ also equipped with the $k$-weighted norm $\left\|\cdot\right\|_{V(D)}$ for all $D\subset\Omega^{\epsilon}$. Throughout this paper, we denote $(\cdot,\cdot)_{D}$ as the $L^{2}(D;\mathbb{C})$ inner product for any Lipschitz domain $D$, $\mathrm{Re}\{\cdot\}$, $\mathrm{Im}\{\cdot\}$ and $\bar{\cdot}$ the real part, the imaginary part and the conjugate of a complex value.

The numerical approach we advance uses the weak formulation for problem (1.1) which is to find $u\in V$ such that

[TABLE]

Here, the sesquilinear form $a:V\times V\to\mathbb{C}$ has the form

[TABLE]

The following properties of the sesquilinear form $a:V\times V\to\mathbb{C}$ play a critical role, these can be found, e.g., in [31, Theorem 3.2 and Corollary 3.3]:

Theorem 2.1** (Properties of the sesquilinear form $a:V\times V\to\mathbb{C}$).**

The following properties hold:

The sesquilinear form $a:V\times V\to\mathbb{C}$ is bounded: There exists a wavenumber $k$ independent constant, ${\rm C}_{\mathrm{b}}$, satisfying:

[TABLE]

• 2.

The following Gårding’s inequality holds:

[TABLE]

The well-posedness of Problem (2.2) can be found, e.g., in [29, Proposition 8.1.3]. Furthermore, there exists some constant, ${\rm C}_{\mathrm{ap}}(k)$, that may depend on the wavenumber $k$ and also on the perforated domain $\Omega^{\epsilon}$, such that the unique solution $u\in V$ to Problem (2.2) fulfills

[TABLE]

Next, we introduce the dual problem to Problem (2.2). For any $w\in L^{2}(\Omega^{\epsilon};\mathbb{C})$, let $z\in V$ be

[TABLE]

then

[TABLE]

2.1 Ansatz space

Since the Gårding’s inequality in Theorem 2.1, combined with the approximation properties of an ansatz space, implies the quasi-optimality of the conforming Galerkin formulation, we now introduce the basic construction of the ansatz space.

Let $\mathcal{T}_{H}$ be a regular partition of the domain $\Omega$ into finite elements with a mesh size $H$. We refer to this partition as coarse grids, and the produced elements as the coarse elements. For each coarse element $K\in\mathcal{T}_{H}$, $K\cap\Omega^{\epsilon}$ is further partitioned into a union of connected fine grid blocks. The fine-grid partition is denoted by $\mathcal{T}_{h}$ with $h$ being its mesh size. Over the fine mesh $\mathcal{T}_{h}$, let $V_{h}$ be the conforming piecewise linear finite element space:

[TABLE]

where $\mathcal{P}_{1}(T)$ denotes the space of linear polynomials on the fine element $T\in\mathcal{T}_{h}$. Then the fine-scale solution $u_{h}\in V_{h}$ satisfies

[TABLE]

The GMsFEM, with which our WEMsFEM shares features and builds from, aims at solving Problem (2.6) on the coarse mesh $\mathcal{T}_{H}$ cheaply, whilst simultaneously maintaining a certain accuracy as compared to the fine-scale solution $u_{h}$. To describe the GMsFEM, we need some notation: The vertices of $\mathcal{T}_{H}$ are denoted by $\{O_{i}\}_{i=1}^{N}$, with $N$ being the total number of coarse nodes. The coarse neighborhood associated with the node $O_{i}$ is denoted by

[TABLE]

The overlap constant $C_{\mathrm{ov}}$ is defined by

[TABLE]

We refer to Figure 2 for an illustration of neighborhoods and elements subordinated to the coarse discretization $\mathcal{T}_{H}$. Throughout, we use $\omega_{i}$ to denote a coarse neighborhood.

Next, we outline the GMsFEM with a conforming Galerkin (CG) formulation. We denote by $\omega_{i}$ the support of the multiscale basis functions. These basis functions are denoted by $\psi_{j}^{\omega_{i}}$ for $j=1,\cdots,\ell_{i}$ for some $\ell_{i}\in\mathbb{N}_{+}$, which is the number of local basis functions associated with $\omega_{i}$. Throughout, the superscript $i$ denotes the $i$-th coarse node or coarse neighborhood $\omega_{i}$. Generally, the GMsFEM utilizes multiple basis functions per coarse neighborhood $\omega_{i}$, and the index $j$ represents the numbering of these basis functions. In turn, the CG multiscale solution $u_{\text{ms}}$ is sought as $u_{\text{ms}}=\sum_{i,j}c_{j}^{i}\psi_{j}^{\omega_{i}}$. Once the basis functions $\psi_{j}^{\omega_{i}}$ are identified, the CG global coupling is given through the variational form

[TABLE]

Here, $V_{\text{ms}}:=\text{span}\{\psi_{j}^{\omega_{i}}:i=1,2,\ldots,N,j=1,2,\ldots,\ell_{i}\}$ denotes the ansatz space.

3 WEMsFEM

We now present our main multiscale method to efficiently solve Problem (1.1) and note that we are especially interested in the cases where the size of the cell $\epsilon$ is very small and the wavenumber $k$ is large. Our method is based on the edge multiscale method proposed in [14]. The main idea is to utilize the wavelets to approximate $u|_{\partial\omega_{i}}$ and then transfer this approximation property into the error over the global domain $\Omega^{\epsilon}$. For the completeness of the presentation, we introduce the wavelets in the following section, the details of which are also found, e.g., in [14].

3.1 Haar wavelet

Let the level parameter and the mesh size be $\ell$ and $h_{\ell}:=2^{-\ell}$ with $\ell\in\mathbb{N}$, respectively, then the grid points on level $\ell$ are

[TABLE]

Let the scaling function $\phi(x)$ and the mother wavelet $\psi(x)$ be given by

[TABLE]

By means of dilation and translation, the mother wavelet $\psi(x)$ can result in an orthogonal decomposition of the space $L^{2}(I)$ with $I:=[0,1]$. To this end, we can define the basis functions on level $\ell\geq 1$ by

[TABLE]

The subspace of level $\ell$ is

[TABLE]

and we note that subspace $W_{\ell}$ is orthogonal to $W_{\ell^{\prime}}$ in $L^{2}(I)$ for any two different levels $\ell\neq\ell^{\prime}$. We denote the subspace in $L^{2}(I)$, up to level $\ell$, by $V_{\ell}^{\mathrm{I}}$ defined by

[TABLE]

The orthogonality of the subspaces $W_{\ell}$ on different levels leads to the relation

[TABLE]

and, consequently, yields the hierarchical structure of the subspace $V_{\ell}$, namely,

[TABLE]

Furthermore, the following orthogonal decomposition of the space $L^{2}(I)$ holds

[TABLE]

Note that one can derive the hierarchical decomposition of the space $L^{2}(I^{d})$ for $d>1$ by means of the tensor product. The following approximation property holds [14, Proposition 3.1]:

Proposition 3.1** (Approximation properties of the hierarchical space $V_{\ell}$ in [14]).**

Let $P_{\ell}$ be $L^{2}(I)$-orthogonal projection onto $V_{\ell}$ for each level $\ell\geq 0$ and let $s>0$. Then there holds

[TABLE]

3.2 The method

We propose our main multiscale algorithm in this section, specifically the corresponding multiscale basis functions are defined locally on each coarse neighborhood independently, and thereby are calculated in parallel. Essentially, we apply wavelets to approximate the solution restricted on each coarse edge. To obtain conforming global basis functions, we utilize the Partition of Unity finite element method [12, 28]; the main idea is to seek local multiscale basis functions in each coarse neighborhood, having certain approximation properties to the exact solution restricted on each coarse neighborhood, and use the fact that the global multiscale basis functions, obtained from those local multiscale basis functions by the partition of unity functions, inherit these approximation properties.

To this end, we begin with an initial coarse space $V^{\text{init}}_{0}=\text{span}\{\chi_{i}\}_{i=1}^{N}$, with the $\chi_{i}$ as the standard multiscale basis functions on each coarse element $K\in\mathcal{T}_{H}$ defined via

[TABLE]

Here $g_{i}$ is affine over $\partial K$ with $g_{i}(O_{j})=\delta_{ij}$ for all $i,j=1,\cdots,N$ and we recall that $\{O_{j}\}_{j=1}^{N}$ are the set of coarse nodes on $\mathcal{T}_{H}$.

Algorithm 1 proceeds as follows: We first construct the local multiscale basis functions on each coarse neighborhood $\omega_{i}$. Given level parameter $\ell\in\mathbb{N}$, and the four coarse edges $\Gamma_{i,k}$ with $k=1,2,3,4$, i.e., $\cup_{k=1}^{4}\Gamma_{i,k}=\partial\omega_{i}$, we let $V_{\ell,k}^{i}$ be the Haar wavelet up to level $\ell$ on the coarse edge $\Gamma_{i,k}$. Introducing $V_{i,\ell}:=\oplus_{k=1}^{4}V_{\ell,k}^{i}$ to be the edge basis functions on $\partial\omega_{i}$, then $V_{i,\ell}$ becomes a good approximation space of dimension $2^{\ell+2}$ to the trace of the solution over $\partial\omega_{i}$, i.e., $u|_{\partial\omega_{i}}$.

Subsequently, we calculate the local multicale basis functions on each coarse neighborhood $\omega_{i}$ with all possible Dirichlet boundary conditions in $V_{i,\ell}$, and denote the resulting local multiscale space as $\mathcal{L}^{-1}_{i}(V_{i,\ell})$ in Step 2. We can then define the global multiscale space as $V_{\text{ms},\ell}^{\text{EW}}$ and obtain the multiscale solution $u_{\text{ms},\ell}^{\text{EW}}$ in Steps 3 and 4.

Note that all the global multiscale basis functions in $V_{\text{ms},\ell}^{\rm EW}$ fulfill the interface condition, i.e. they satisfy the Neumman boundary condition on every interface of every inclusion, as in the second equation in (1.1). This is due to the construction of the partition of unity functions $\chi_{i}$ in (3.1).

3.3 Local projection

An important element of our algorithm is the knowledge of the approximation properties to the edge basis functions $V_{i,\ell}:=\oplus_{k=1}^{4}V_{\ell,k}$ on each coarse neighborhood $\partial\omega_{i}$.

Let the $L^{2}(\partial\omega_{i})$-orthogonal projection $\mathcal{P}_{i,\ell}$ onto the local multiscale space up to level $\ell$: $L^{2}(\partial\omega_{i})\to V_{i,\ell}$ be defined by

[TABLE]

Here, we denote $\psi_{j}$ for $j=1,\cdots 2^{\ell+2}$ as the Haar wavelet defined on the four edges of $\omega_{i}$ of level $\ell$ and the local operator $\mathcal{L}_{i}$ is defined in Algorithm 1.

Let $\text{diam}(\Omega^{\epsilon})$ be the diameter of the bounded domain $\Omega^{\epsilon}$. Define

[TABLE]

Then the positive constants ${\rm C}_{\mathrm{poin}}(\omega_{i})$ and ${\rm C}_{\mathrm{poin}}(\Omega^{\epsilon})$ are independent of the wavenumber $k$ and the coarse mesh $\mathcal{T}_{H}$. Note that we will utilize the same constant ${\rm C}_{\mathrm{poin}}(\omega_{i})$ to denote the constant from the Poincaré inequality.

Assumption 3.1** (Scale Resolution Assumption).**

We assume that the coarse mesh size $H$ is sufficiently small to satisfy the following inequality

[TABLE]

To simplify the notation, we denote

[TABLE]

Remark 3.1**.**

Similar resolution assumption as Assumption 3.1 is commonly seen in the literature, e.g., [31]. If we take the coarse scale mesh grid $H\ll k^{-1}$, then ${\rm C}_{\mathrm{est}}\approx 1$.

4 Convergence rate of Algorithm 1

The convergence rate of this algorithm is clearly an important detail and, perhaps remarkably, this can be obtained. The proof is inspired by the techniques developed in [14, 21], where a local decomposition of the solution $u$ restricted on each coarse neighborhood $\omega_{i}$ for $i=1,2,\ldots,N$, namely, $u|_{\omega_{i}}$, is introduced. We also analyze the local approximation properties of the multiscale basis functions constructed in Step 2, Algorithm 1 to each component of the decomposition of $u|_{\omega_{i}}$ in Section 4.1. Subsequently, the global approximation properties of the ansatz space $V_{\text{ms},\ell}^{\rm EW}$ and the convergence of Algorithm 1 are investigated in Section 4.2.

4.1 Local decomposition of the solution

The solution $u$ satisfies the following equation

[TABLE]

which can be split into three parts, namely

[TABLE]

Here, the three components $u^{i,\mathrm{I}}$, $u^{i,\mathrm{II}}$ and $u^{i,\mathrm{III}}$ are respectively given by

[TABLE]

[TABLE]

and

[TABLE]

Here, \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\omega_{i}}v:=|\omega_{i}|^{-1}\int_{\omega_{i}}v\,\mathrm{d}x denotes the average of the function $v\in L^{1}(\omega_{i})$ over each coarse neighborhood $\omega_{i}$. Recall that $v^{i}$ is defined in Step 2 of Algorithm 1.

We first show that $u^{i,\mathrm{I}}$ is negligible thanks to the local basis function $v^{i}$:

Lemma 4.1**.**

Let the Scale Resolution Assumption 3.1 be valid. Let $u^{i,\mathrm{I}}$ be defined in (4.2) and let $f\in L^{2}(\Omega^{\epsilon})$. Then there holds

[TABLE]

Proof.

Multiplying (4.2) by $u^{i,\mathrm{I}}$, and taking its integral over $\omega_{i}$, leads to

[TABLE]

wherein an application of Hölder’s inequality and the Poincaré inequality proves

[TABLE]

After moving the first term on the right of the previous estimate to the left, and noting the Scale Resolution Assumption 3.1, we obtain

[TABLE]

Finally, an application of the Poincaré inequality and the Scale Resolution Assumption 3.1 shows the desired assertion. ∎

Recalling that $z\in V$ is the solution to Problem (2.4) for given $w\in L^{2}(\Omega^{\epsilon};\mathbb{C})$, analogously to Decomposition (4.1), the following decomposition is valid:

[TABLE]

Then a similar argument as in Lemma 4.1 leads to the following estimate:

Lemma 4.2**.**

Let the Scale Resolution Assumption 3.1 hold. Let $z^{i,\mathrm{I}}$ be defined in (4.3) and let $w\in L^{2}(\Omega^{\epsilon})$. Then it holds that

[TABLE]

Remark 4.1** (On the decomposition (4.1) and (4.3)).**

Following from Lemmas 4.1 and 4.2, only one multiscale interial basis function in Step 2 of Algorithm 1 is sufficient if $\mathcal{O}(H)$ convergence rate is desired or sufficient. Otherwise, one can construct extra local multiscale basis functions to approximate the first component $u^{i,\mathrm{I}}$.

Since $u^{i,\mathrm{III}}$ is of rank-one, which can be represented with one multiscale basis function $v^{i}$ in Decomposition (4.1). Lemma 4.1 indicates that we need only construct a proper ansatz space for the second part $u^{i,\mathrm{II}}$ to ensure a good ansatz space for $u|_{\omega_{i}}$. We now prove that the multiscale basis functions, constructed in Step 2 of Algorithm 1, span an appropriate ansatz space with good approximation properties:

Theorem 4.1** (Approximation properties of the projection $\mathcal{P}_{i,\ell}$).**

Let Assumption 3.1 hold and let $e\in V(\Omega^{\epsilon}\cap\omega_{i})$ satisfy

[TABLE]

Then it holds that

[TABLE]

Proof.

We can obtain from Theorem A.1:

[TABLE]

Whereas an application of [14, Eqn. (5.6)] leads to

[TABLE]

We only need to estimate $\left|u^{i,\mathrm{II}}\right|_{H^{1}\left(\omega_{i}\right)}$. Note that $r:=u-u^{i,\mathrm{II}}$ satisfies

[TABLE]

A similar estimate as used in the proof to Lemma 4.1 shows that

[TABLE]

Finally, an application of the triangle inequality proves

[TABLE]

and this, together with (4.5), proves the desired assertion. ∎

4.2 Approximation properties of the ansatz space $V_{\text{ms},\ell}^{\rm EW}$

The approximation properties of the ansatz space $V_{\text{ms},\ell}^{\rm EW}$ follow from the theorem below:

Theorem 4.2** (Approximation properties of the multiscale space $V_{\text{ms},\ell}^{\rm EW}$).**

Let Assumption 3.1 hold and assume that $f\in L^{2}(\Omega^{\epsilon})$. Let $\ell\in\mathbb{N}_{+}$ and $u\in V$ be the solution to Problem (1.1). There then holds

[TABLE]

Furthermore, let $w\in L^{2}(\Omega^{\epsilon})$ and let $z\in V$ be the solution to Problem (2.4). Then it holds

[TABLE]

Here, ${\rm C}_{\mathrm{\ref{thm:last}}}$ is a positive constant independent of the wavenumber $k$ or coarse mesh size $H$.

Proof.

We will only prove the first assertion (4.6) since the second assertion (4.7) can be derived in a similar manner.

Recall the local decomposition in (4.1) on each coarse neighborhood $\omega_{i}$ for $i=1,\ldots,N$. Let

[TABLE]

We prove that $v$ is a good approximation to $u$.

Denote ${e}_{\mathrm{}}:=u-v$. Then the property of the partition of unity of $\{\chi_{i}\}_{i=1}^{N}$ leads to

[TABLE]

Taking its squared energy norm, and using the overlap condition (2.8), we arrive at

[TABLE]

Then Young’s inequality implies

[TABLE]

Using the product rule, and the Poincaré inequality, we obtain

[TABLE]

Then Lemma 4.1 yields

[TABLE]

Analogously, we can derive the following upper bound for the second term by Theorem 4.1:

[TABLE]

Inserting these two estimates into (4.2), and utilizing the overlapping condition (2.8), leads to

[TABLE]

Furthermore, inserting (2.3) into (4.10) shows (4.6). This completes the proof. ∎

Finally, we are ready to present the main result of this section:

Proposition 4.1** (Error estimate for Algorithm 1).**

Assume that $f\in L^{2}(\Omega^{\epsilon})$ and let the coarse mesh size $H$ and the level parameter $\ell\in\mathbb{N}_{+}$ satisfy

[TABLE]

Let $u\in V$ and $u_{H}\in V_{\rm{ms},\ell}^{\rm EW}$ be the solution to Problem (1.1), and from Algorithm 1, respectively. There holds

[TABLE]

Proof.

Since $V_{\text{ms},\ell}^{\rm EW}\subset V$, Gårding’s inequality in Theorem 2.1 implies

[TABLE]

for all $v_{\text{ms},\ell}^{\text{EW}}\in V_{\text{ms},\ell}^{\text{EW}}$.

Next we estimate $\left\|u-u_{\text{ms},\ell}^{\text{EW}}\right\|_{L^{2}\left(\Omega^{\epsilon}\right)}$: Let $z\in V$ be the solution to

[TABLE]

then for all $z_{\text{ms},\ell}^{\text{EW}}\in V_{\text{ms},\ell}^{\text{EW}}$ we obtain that

[TABLE]

Furthermore, an application of Theorem 2.1 leads to

[TABLE]

By (4.7), and condition (4.11), we arrive at

[TABLE]

This, together with (4.13), leads to

[TABLE]

Consequently, we obtain

[TABLE]

Finally, an application of the boundedness of the sesquilinear form $a(\cdot,\cdot)$ Theorem 2.1 and the approximation property (4.12) prove the desired assertion. ∎

5 Numerical experiments

We present numerical experiments to show the performance of Algorithm 1 with locally periodic perforations in Sections 5.1 and 5.2, and with random perforations in Section 5.3. We consider two different source excitations placed at $(0,0)$ and $(1.44,0)$ for each model, and take the wavenumber $k=64$.

The coarse mesh $\mathcal{T}_{H}$ is a regular partition of the domain $\Omega^{\epsilon}$ into finite elements with a mesh size $H$. Then each coarse element is further partitioned into a union of connected fine grid blocks. The fine-grid partition is denoted by $\mathcal{T}_{h}$, which provides a sufficiently fine mesh for standard finite element solvers to get a reference solution; for sufficient accuracy we take $h:={1}/{320}$. In addition, the Perfectly Matched Layer (PML) is utilized to absorb the outgoing wave, see Appendix B for more details.

As is usual in the finite element literature we use the $L^{2}(\Omega^{\epsilon})$ and $H^{1}(\Omega^{\epsilon})$-relative errors to assess the accuracy of the scheme, which are defined by

[TABLE]

5.1 Performance of Algorithm 1: level

parameter $\ell$

We consider in this section the two perforated models as shown in Figure 1. To describe the computational domain $\Omega^{\epsilon}$, we use Eq.(2.1). Let $\Omega:=[-2.4,2.4]^{2}$, $\Omega_{1}:=[-1,1]^{2}$ and the size of the cell $\epsilon:=1/6$. The perforations in a unit cell are $Q:=[0.25,0.75]^{2}$ and $Q:=[0.1,0.9]^{2}$ in models 1 and 2. Note that there are $144$ perforations in both models, this strong heterogeneity in the computational domain $\Omega^{\epsilon}$ makes Problem (1.1) even harder.

Figure 3 demonstrates dynamic anisotropy, also known as self-collimation, [3, 34] which is a striking, frequency sensitive and dependent, effect. Naively, one might assume that wave excitation of the crystal, at its centre, would lead to isotropic wavefronts within the crystal; this is indeed the case in the standard homogenisation, low-frequency long-wave, limit where the wavelength is much larger than the cell-to-cell spacing. However, as the frequency increases and wavelength decreases, Bragg scattering occurs with constructive or destructive interference leading to well-defined frequency windows (band-gaps) within which wave propagation is disallowed. Interference also occurs that acts to create anisotropic wavefronts with the most severe example being that where all the wave energy is channeled in specific directions; in terms of homogenisation there are variants that are developed for high-frequencies [7] that show the effective medium PDE changing its character from elliptic to hyperbolic with these directions of self-collimation being the characteristics of the hyperbolic system [26]. Figure 4 shows lensing, another effect created by the crystal whereby a partial image of the source (on the right) forms to the left of the crystal.

As we can see from Figures 3 and 4, even with $\ell=2$, one can observe the wave scattering phenomenon resulting from the perforated structure with clear agreement between the multiscale solution and the reference solution with full capture of the microscale feature of the wavefield; we now quantify this agreement. We test the convergence of Algorithm 1 with respect to the level parameter $\ell$ for both of the perforated domains shown in Figure 1. To this end, we fix the coarse scale mesh size $H:=1/10$. Recall that the level parameter $\ell$ determines the number of multiscale basis functions in each coarse neighborhood $\omega_{i}$ for $i=1,\cdots,N$ with $N$ as the number of coarse grids; specifically this number is $2^{\ell+2}+1$ and the level parameter $\ell$ shows the complexity of Algorithm 1.

Figure 5 shows the $L^{2}(\Omega^{\epsilon})$ and $H^{1}(\Omega^{\epsilon})$-relative errors versus $\ell$, and both the $L^{2}(\Omega^{\epsilon})$ and $H^{1}(\Omega^{\epsilon})$-relative errors decay rapidly as more wavelet basis functions are added. For example, for the case that the source lies at the center shown in Figure 5(a), the $L^{2}(\Omega^{\epsilon})$ errors decrease from $120.0\%$ to $4.08\%$ as the level parameter $\ell$ increases from [math] to $3$. Figure 5 suggests that Algorithm 1 with level parameter $\ell=2$ yields an accurate solver with little gain from going to higher $\ell$.

Similar performance for the perforated domain $\Omega^{\epsilon}$ depicted in model 2 is obtained, see Figures 6 and 7. The second model as compared with the perforated domain in model 1 has the gap between perforations much narrower, we expect stronger singularities and finer structure in the solution, and this provides insight on method robustness. In terms of the physics, the frequency has remained fixed and altering the perforation size alters the dynamic anisotropy slightly, we observe strong directionality and concentration of the highly oscillatory wave field in the narrow gaps.

The relative error decay history is shown in Figure 8 and we find similar relative error decay behavior as in Figure 5 and Algorithm 1 is both efficient and accurate. For instance, for the case that the source lies in the center, the $L^{2}(\Omega^{\epsilon})$-relative error reaches below $10\%$ even with the level parameter $\ell=1$.

5.2 Performance of Algorithm 1:

coarse-scale mesh size $H$

Earlier we established theoretically that the error induced by Algorithm 1 can attain $\mathcal{O}(H)$ in Proposition 4.1, upon the condition on the coarse mesh size $H$ and the level parameter $\ell$, cf. (4.11). We now test how the algorithm performs with respect to a different, finer, coarse-scale mesh, we take its size $H:=1/20$, so that the coarse mesh $\mathcal{T}_{H}$ will cross the perforations.

The error decay history for the two different perforated domains, depicted in Figure 1 with centered and right sources, are plotted in Figures 9 and 10, respectively. There is error decay as the level parameter $\ell$ with very rapid decay and we conclude that multiscale solutions with sufficient accuracy are achieved with the level parameter $\ell=2$. Comparing with the mesh of $H:=1/10$, Figures 5 and 8, we conclude that as the coarse-scale mesh size $H$ decreases, the performance of Algorithm 1 significantly improves; this agrees with the predictions of Proposition 4.1.

5.3 Performance of Algorithm 1: random

perforations

The motivation for the development of the multiscale WEMsFEM was to emerging problems in finite periodic crystals, but the methodology is not reliant on any periodicity and we investigate the algorithm’s performance in more general situations. A natural case to consider is that of random perforated domains as in Figure 11. The size of the perforation for model 3 is 0.08, and for model 4 is 0.24 and we consider $H=1/20\text{ and }1/10$ for model 3 and $H=1/10$ for model 4.

Firstly, we present the reference solution and multiscale solution from Algorithm 1 with model 3 as the perforated domain, a centered source term, a coarse mesh size $H:=1/10$ and level parameter $\ell=2$ in Figure 12. The $L^{2}(\Omega^{\epsilon})$-relative error is 3.70%; further decreasing the coarse mesh size $H$ or increasing the level parameter $\ell$ improves the accuracy as shown in Figure 13.

Next, we study the performance of Algorithm 1 in the perforated domain of model 4. The perforations in model 4 cross the neighborhood boundary and we depict the resulting four coarse neighborhood $\omega_{i}$ in Figure 14. The local solvers in Algorithm 1 are now defined in some L-shaped domains and consequently, the perforated domain of model 4 is much more challenging numerically.

Nonetheless, we observe similar performance as for model 3, with a centered source, $H:=1/10$ and $\ell=2$ we show the results in in Figure 15; the corresponding $L^{2}(\Omega^{\epsilon})$-relative error is 3.82%. Analogously to model 3, further increases in the level parameter $\ell$, or decrease in the coarse grid size $H$, further improve the performance of our algorithm, see Figure 16 for more details.

6 Conclusion

We demonstrate that the Wavelet-based Edge Multiscale Finite Element Method (WEMsFEM) for Helmholtz problems in perforated domains, with possibly large wavenumbers, are an effective alternative to standard Finite Element methods with advantages for multiscale problems. Such problems have many applications in photonic and phononic crystals and having efficient, faster, alternatives to the existing finite element approaches is much needed particularly for crystals created from many cells and operating at high frequencies.

We have created both the required theory, and error estimates, and then tested the convergence analysis and numerical performance of the algorithm against examples of physical interest. Under the usual resolution assumption that the product of the coarse-scale mesh size $H$ and the wavenumber $k$ is bounded above by a certain constant and the level parameter $\ell$ is sufficiently large, we prove $\mathcal{O}(H)$ convergence of our methods. Our theoretical results are supported by extensive 2-d numerical simulations. The success of this two-dimensional study has motivated further practical tests of this algorithm for 3-d Helmholtz problems in perforated domains and these are currently under investigation.

Appendix A Very-weak solutions to the Helmholtz problem

We establish in this section an a priori estimate, cf. Theorem A.1, which is utilized in the proof to Theorem 4.1. Throughout this section, $\omega_{i}$ is one coarse neighborhood as defined in (2.7) for all $i=1,2,\cdots,N$.

Let $g\in L^{2}(\Omega^{\epsilon}\cap\omega_{i})$, and $v\in H^{1/2}(\Omega^{\epsilon}\cap\omega_{i})$ be the solution to the following problem:

[TABLE]

[TABLE]

This test space $X(\omega_{i})$ is endowed with the norm $\|\cdot\|_{X(\omega_{i})}$:

[TABLE]

Then we propose the following weak formulation corresponding to Problem (A.1): seeking $v\in L^{2}(\Omega^{\epsilon}\cap\omega_{i})$ such that

[TABLE]

Theorem A.1**.**

Let the Scale Resolution Assumption 3.1 be valid. Given $g\in L^{2}(\Omega^{\epsilon}\cap\omega_{i})$. Let $v$ be the solution to (A.1). Then there holds

[TABLE]

To prove Theorem A.1, one has to first derive the $L^{2}(\partial\omega_{i}\backslash\partial Q_{1}^{\epsilon})$-estimate of the normal trace $\frac{\partial z}{\partial n}$ for any $z\in X(\omega_{i})$. This is established in the following theorem:

Theorem A.2**.**

Let the Resolution Assumption 3.1 be valid. Let $w\in L^{2}(\Omega^{\epsilon}\cap\omega_{i})$ and let $z\in X(\omega_{i})$ satisfy

[TABLE]

Then there holds

[TABLE]

Here, ${\rm C}_{\mathrm{weak}}$ is a positive constant independent of the wavenumber $k$ and the mesh size $H$, which can change values among equations.

Proof.

An application of the Poincaré inequality implies

[TABLE]

Testing (A.4) with $z$ and applying the Poincaré inequality, together with the former result, we obtain

[TABLE]

Thanks to the Resolution Assumption (3.1), this yields

[TABLE]

Furthermore, by (A.5), we arrive at

[TABLE]

One the other hand, a direct calculation results in

[TABLE]

This, together with (A.7) and the Resolution Assumption (3.1), yields

[TABLE]

Note that the interface $\partial Q_{1}^{\epsilon}\cap\omega_{i}$ has sufficient smoothness, we have the following a priori estimate

[TABLE]

This, together with (A.6) and applying interpolation between $H^{1}(\Omega^{\epsilon}\cap\omega_{i})$ and $H^{2}(\Omega^{\epsilon}\cap\omega_{i})$ yields the $H^{3/2}(\Omega^{\epsilon}\cap\omega_{i})$ regularity estimate

[TABLE]

Since differentiation is continuous from $H^{3/2}(\Omega^{\epsilon}\cap\omega_{i})$ to $H^{1/2}(\Omega^{\epsilon}\cap\omega_{i})$, by the trace theorem, we have

[TABLE]

which, together with (A.8) and the Resolution Assumption 3.1, proves the desired assertion. ∎

Proof to Theorem A.1.

The first result can be proved in a similar manner as [21, Theorem A.1], with the help of Theorem A.2.

To prove the second assertion, recall that $\chi_{i}$ is the bilinear function supported in $\omega_{i}$ and $\chi_{i}=0$ on $\partial\omega_{i}\backslash\partial Q_{1}^{\epsilon}$. Multiplying (A.1) by $\chi_{i}^{2}v$ and applying integration by parts, we arrive at

[TABLE]

Then an application of the Young’s inequality implies

[TABLE]

After taking the square root over the previous estimate, utilizing the first assertion and the Scale Resolution Assumption 3.1, the second assertion is proved. ∎

Appendix B PML for the Helmholtz equation

To effectively absorb the outgoing wave, we adopt the Perfectly Matched Layer (PML) [1, 4] in our implementation. Without loss of generality, let the computation domain including the PML areas be $D=(0,1)^{2}$. We follow the notations in [13], we define

[TABLE]

and

[TABLE]

and

[TABLE]

where $x_{1}$ and $x_{2}$ are the space variables. $\xi$ is the thickness of the PML. The PML method is to replace $\partial_{1}$ with $g_{1}(x_{1})\partial_{1}$ and $\partial_{2}$ with $g_{2}(x_{2})\partial_{2}$, respectively. Then then Equation (1.1) becomes

[TABLE]

In our implementation, we take $C:=100$ and the thickness of the PML $\xi$ equals to one wavelength.