TY - JOUR
AU1 - Li, Guanglian
AU2 - Peterseim, Daniel
AU3 - Schedensack, Mira
AB - Abstract We formulate a stabilized quasi-optimal Petrov–Galerkin method for singularly perturbed convection–diffusion problems based on the variational multiscale method. The stabilization is of Petrov–Galerkin type with a standard finite element trial space and a problem-dependent test space based on pre-computed fine-scale correctors. The exponential decay of these correctors and their localization to local patch problems, which depend on the direction of the velocity field and the singular perturbation parameter, are rigorously justified. Under moderate assumptions, this stabilization guarantees stability and a quasi-optimal rate of convergence for arbitrary mesh Péclet numbers on fairly coarse meshes at the cost of additional inter-element communication. 1. Introduction Given a domain $${\it{\Omega}}\subset\mathbb{R}^2$$, a singular perturbation parameter $$0<\epsilon\leq 1$$, a velocity field $$b\in (L^\infty({\it{\Omega}}))^2$$ and some force $$f\in H^{-1}({\it{\Omega}})$$, the convection–diffusion equation seeks $$u\in V:=H^1_0({\it{\Omega}})$$ such that −ϵΔu+b⋅∇u =fin Ω,u =0on ∂Ω. (1.1) We assume that the velocity field $$b$$ is incompressible, i.e., $$\nabla\cdot b=0$$. The focus of this article is the convection-dominated regime with large Péclet number $$\mathrm{Pe}=\left\|{b}\right\|_{L^{\infty}\left({{\it{\Omega}}}\right)}/\epsilon$$. For reasonably small Péclet numbers, classical Galerkin finite element methods (FEMs) perform well. However, if the Péclet number increases then steep gradients of $$u$$ occur and boundary layers appear, requiring a much finer mesh to capture the characteristic width of these boundary layers. Consequently, local corrections are needed at the layers and a numerical method in which the smooth solution regions are not polluted by the layers is desirable. The thickness of the parabolic layer is $$\mathcal{O}(\sqrt{\epsilon})$$ and $$\mathcal{O}(\epsilon)$$ for the exponential layer, which have to be resolved for a stable approximation with a standard Galerkin FEM. Furthermore, it holds that $$\left|{u}\right|_{H^{1}\left({{\it{\Omega}}^*}\right)}=\mathcal{O}(\epsilon^{-{1}/{4}})$$ and $$\left|{u}\right|_{H^{1}\left({{\it{\Omega}}^o}\right)}=\mathcal{O}(\epsilon^{-{1}/{2}})$$ with small neighborhoods $${\it{\Omega}}^*$$ and $${\it{\Omega}}^o$$ of the parabolic and the exponential boundary layer, respectively (Roos et al., 2008; John & Schmeyer, 2009). Numerous numerical methods have been proposed in the past few decades aimed at solving the convection-dominated problem (1.1) efficiently and accurately. Upwinding methods for stabilization of the exponential boundary layers combined with refinement near the parabolic boundary layers are formulated. Among them are the streamline upwind/Petrov–Galerkin (SUPG) method or Galerkin least squares (GLS) method (Franca et al., 1992; Christiansen et al., 2016), $$hp$$ finite element methods (Melenk, 1997, 2002), discontinuous Petrov–Galerkin methods (Demkowicz et al., 2012), residual-free bubble (RFB) approaches (Brezzi et al., 2000; Cangiani & Süli, 2005a,b), methods with an additional nonlinear diffusion (Barrenechea et al., 2016), methods with stabilization by local orthogonal sub scales (Codina, 2000) and hybridizable discontinuous Galerkin methods (Qiu & Shi, 2015). Among the multiscale methods are variational multiscale (VMS) methods (Hughes & Sangalli, 2007; Larson & Målqvist, 2009), multiscale finite element methods (Park & Hou, 2004; Calo et al., 2016), multiscale hybrid-mixed methods (Harder et al., 2015) and local orthogonal decomposition (LOD) methods (Elfverson, 2015). Specifically, the residual-based stabilization methods (SUPG, GLS and RFB) incorporate global stability properties into high accuracy in local regions away from the boundary layers. We refer to Roos et al. (2008) for an overview of robust numerical methods for singular perturbed problems. In this article, our focus is on the construction and the error analysis of a stable and accurate LOD method based on Hughes & Sangalli (2007), Peterseim (2017) and Målqvist & Peterseim (2011). VMS methods were designed for solving multiscale problems by embedding fine-scale information into a coarse-scale framework. Essentially, the efficiency and accuracy rely on the construction of a problem-dependent stable projector from a larger fine space onto a relatively much smaller coarse space. Our motivation for this article originates from Hughes & Sangalli (2007), where the authors derived an explicit formula for the one-dimensional fine-scale Green’s function arising in VMS methods. The smaller the support of the fine-scale Green’s function, the more favorable the localized method (e.g., Målqvist & Peterseim, 2011; Peterseim, 2017) in solving (1.1). In particular, the authors compared the fine-scale Green’s functions derived by the $$L^2$$-projector with that derived by the $$H_{0}^{1}$$-projector and concluded that the latter outweighed the former in the one-dimensional case. In addition, examples were shown for the two-dimensional case that the $$H_{0}^{1}$$-projector would exceed the $$L^2$$-projector as well. There is a recent work (Elfverson, 2015) on the convection–diffusion problem employing the $$L^2$$-projector in the framework of the VMS and LOD methods. The author shows convergence of the localized method and tests the method using the $$H_{0}^1$$-projector and claims that the superiority of the $$H_{0}^{1}$$-projector over the $$L^2$$-projector is not valid for the two-dimensional case. This claim is supported by Fig. 2 below. Fig. 1. View largeDownload slide Standard nodal basis function $$\lambda_z$$ with respect to the coarse mesh $${\mathcal{T}}_H$$ (top left), corresponding ideal corrector $${\mathcal{C}} \lambda_z$$ (top right) and corresponding test basis function $$(1-{\mathcal{C}} )\lambda_z$$ (bottom left). The bottom right figure shows a top view of the modulus of the test basis function $$(1-{\mathcal{C}}) \lambda_z$$ with logarithmic color scale to illustrate the exponential decay property. The underlying data are $$b=[\cos(0.7),\sin(0.7)]$$ and $$\epsilon=2^{-7}$$. Fig. 1. View largeDownload slide Standard nodal basis function $$\lambda_z$$ with respect to the coarse mesh $${\mathcal{T}}_H$$ (top left), corresponding ideal corrector $${\mathcal{C}} \lambda_z$$ (top right) and corresponding test basis function $$(1-{\mathcal{C}} )\lambda_z$$ (bottom left). The bottom right figure shows a top view of the modulus of the test basis function $$(1-{\mathcal{C}}) \lambda_z$$ with logarithmic color scale to illustrate the exponential decay property. The underlying data are $$b=[\cos(0.7),\sin(0.7)]$$ and $$\epsilon=2^{-7}$$. Fig. 2. View largeDownload slide Impact of the choice of interpolation operator on the decay of the ideal test basis function $$(1-{\mathcal{C}} )\lambda_z$$ in the under-resolved regime $$H\gg \epsilon$$: nodal interpolation (left), $$L^2$$-projection (middle) and $$H^1_0$$-projection (right). The underlying data are $$b=[\cos(0.7),\sin(0.7)]$$ and $$\epsilon=2^{-7}$$. Fig. 2. View largeDownload slide Impact of the choice of interpolation operator on the decay of the ideal test basis function $$(1-{\mathcal{C}} )\lambda_z$$ in the under-resolved regime $$H\gg \epsilon$$: nodal interpolation (left), $$L^2$$-projection (middle) and $$H^1_0$$-projection (right). The underlying data are $$b=[\cos(0.7),\sin(0.7)]$$ and $$\epsilon=2^{-7}$$. In the one-dimensional case, the $$H^1_0$$-projection equals the nodal interpolation. Therefore, another possible generalization of the one-dimensional case to higher dimensions is to use nodal interpolation in the VMS method. Doing this, we take into account that nodal interpolation is not well defined and stable as an operator on $$H^1$$ functions. After regularization using a very fine reference discretization of scale $$h\lesssim\epsilon$$, the use of nodal interpolation can be justified as the stability constant degenerates with $$\log(H/h)$$ in two dimensions only (see (3.1) below). This approach was previously utilized in Larson & Målqvist (2009) and is shown to work better than averaging-type operators in the regime of large mesh Péclet numbers (see Fig. 2 below). In this article, we show that a VMS method based on the nodal interpolation operator coupled with a Petrov–Galerkin method is stable and locally quasi-optimal for the convection-dominated problem (1.1) with no spurious oscillations and no smearing. As for other elliptic partial differential equations (PDEs), the ideal VMS method is turned into a practical method by localizing the support of the VMS basis functions (Peterseim 2016b). Inspired by the numerical results of the fine-scale Green’s functions displayed in Hughes & Sangalli (2007) and the proof in our article as well, a $$b$$-biased local region is proposed as the numerical domain for approximating the ideal method. The convergence of this localization is proved under the assumption that the local region is sufficiently large. In three dimensions, the method is still applicable, but the instability of nodal interpolation is more severe and the size of the subdomains may be prohibitively large to compensate this negative effect. The remainder of this article is organized as follows. In Section 2, a detailed description of the problem considered in this article is shown. In Section 3, we propose a new VMS method based on the nodal interpolation and denote it as the ideal method. Its stability and local quasi-optimality are displayed. In Section 4, we estimate the error of the global correctors outside a certain local patch and show an exponential decay of the error with respect to the size of the local patch. Inspired by the results in Section 4, we formulate the localization algorithm in Section 5 for the ideal method proposed in Section 3 and display the stability of this algorithm as well as the convergence. A numerical experiment is provided in Section 6 for the validation of our method and we conclude this article with conclusions in Section 7. 2. Model problem and standard finite elements We assume that the parameter $$\epsilon\leq 1$$ and $$b\in L^{\infty}({\it{\Omega}};\mathbb{R}^2)$$ is a divergence-free vector field and we define the bilinear form $$a$$ on $$V\times V$$ associated to (1.1) by a(u,v)=ϵ∫Ω∇u⋅∇vdx+∫Ω(b⋅∇u)vdxfor all u,v∈V. (2.1) Because $$\nabla\cdot b=0$$, an integration by parts implies that the bilinear form $$a$$ is $$V$$-elliptic, i.e., a(v,v) =ϵ|v|H1(Ω)2 for all v∈V. (2.2) Furthermore, a Poincaré–Friedrichs inequality leads to the existence of some constant $$C({\it{\Omega}},b)$$ that may depend on (the diameter of) the domain $${\it{\Omega}}$$ and the $$L^\infty$$ norm of $$b$$ such that $$a$$ is continuous, i.e., for all $$u,v\in V$$ it holds that a(u,v) ≤ϵ|u|H1(Ω)|v|H1(Ω)+‖b‖L∞(Ω)|u|H1(Ω)‖v‖L2(Ω) ≤C(Ω,b)|u|H1(Ω)|v|H1(Ω). (2.3) Here, we used that $$\epsilon\leq 1$$. Throughout this article, $$A\lesssim B$$ abbreviates that there exists a constant $$C>0$$ independent of $$\epsilon$$, $$h$$ and $$H$$ ($$h$$ and $$H$$ will be defined later), such that $$A\leq C B$$, and let $$A\gtrsim B$$ be defined as $$B\lesssim A$$, and $$A\approx B$$ abbreviates $$A\lesssim B\lesssim A$$. We assume that $$b\neq (0\ \ 0)$$$$\left\|{b}\right\|_{L^{\infty}\left({{\it{\Omega}}}\right)}\approx 1$$. Let $$\langle\bullet,\bullet\rangle_{H^{-1}({\it{\Omega}})\times H^1_0({\it{\Omega}})}$$ denote the dual pairing of $$H^{-1}({\it{\Omega}})$$ and $$H^1_0({\it{\Omega}})$$. We consider the variational form of (1.1): {find u∈V such that for all v∈Va(u,v)=⟨f,v⟩H−1(Ω)×H01(Ω). (2.4) By virtue of the $$V$$-ellipticity and $$V$$-continuity of $$a$$ from (2.2) and (2.3) and the Lax–Milgram lemma, problem (2.4) has a unique solution in $$V$$. Let $${\mathcal{T}}_h$$ be a shape-regular triangulation of the domain $${\it{\Omega}}$$, where $$h$$ represents the minimal diameter of all triangles in $${\mathcal{T}}_h$$. Given a triangulation $${\mathcal{T}}$$, let P1(T) :={v∈C0(Ω)∣v|K∈P1(K) for all K∈T} denote the space of piecewise linear finite elements and define $$V_h:={\mathcal{P}}_1({\mathcal{T}}_h)\cap V$$. Let $$u_h\in V_h$$ denote the reference solution, which is defined as the Galerkin approximation that satisfies a(uh,vh)=⟨f,vh⟩H−1(Ω)×H01(Ω)for all vh∈Vh. (2.5) Taking advantage of the ellipticity and continuity of $$a$$ from (2.2) and (2.3) on $$V\times V\supset V_h\times V_h$$, the Lax–Milgram lemma implies that the fine-scale solution $$u_h$$ of (2.5) exists and is unique on $$V_h$$. We assume that $$\epsilon\ll 1$$ is a small parameter and that $${\mathcal{T}}_h$$ resolves $$\epsilon$$ in the sense that $$u_h$$ is a good approximation of $$u$$, e.g., if hmax‖b‖L∞(Ω)/ϵ≲1 (2.6) with the maximal mesh size $$h_\mathrm{max}$$ of $${\mathcal{T}}_h$$. It holds that |u−uh|H1(Ω)≲(1+hmax‖b‖L∞(Ω)ϵ)infvh∈Vh|u−vh|H1(Ω). If, in addition, the solution $$u$$ of (2.4) satisfies $$u\in H^2({\it{\Omega}})$$, standard interpolation estimates lead to |u−uh|H1(Ω)≲hmax(1+hmax‖b‖L∞(Ω)ϵ)‖u‖H2(Ω), with a hidden constant independent of $$\epsilon$$. Note, however, that $$\|u\|_{H^2({\it{\Omega}})}$$ depends on $$\epsilon$$. 3. The ideal method In this section, we introduce a variational multiscale method based on nodal interpolation, which yields a locally best approximation of the reference solution $$u_h\in V_h$$ from (2.5) and which is computed on a feasible coarse underlying mesh $${\mathcal{T}}_H$$. We assume that $${\mathcal{T}}_H$$ is a regular quasi-uniform triangulation of the domain $${\it{\Omega}}$$ with maximal mesh size $$H$$, such that $${\mathcal{T}}_h$$ is a refinement of $${\mathcal{T}}_H$$. Let $$\mathcal{N}_H$$ denote the nodes in $${\mathcal{T}}_H$$ and $$\mathrm{mid}_K$$ the barycenter for each coarse element $$K\in {\mathcal{T}}_H$$. The maximal mesh size $$H$$ of $${\mathcal{T}}_H$$ represents a computationally feasible scale that is typically much larger than $$\epsilon$$. Altogether, the target regime is then 00$$, such that maxK∈TH#{T∈TH∣K⊂ΩT,ℓ,b}≤Col,ℓ(ϵ). (4.4) Fig. 3. View largeDownload slide Element patches $${\it{\Omega}}_{T,\ell,b}$$ for $$b=[\cos(0.7),\sin(0.7)]$$, $$\epsilon=2^{-7}$$ and $$\ell=1,2,3$$ (from left to right) as they are used in the localized corrector problem (5.1). Fig. 3. View largeDownload slide Element patches $${\it{\Omega}}_{T,\ell,b}$$ for $$b=[\cos(0.7),\sin(0.7)]$$, $$\epsilon=2^{-7}$$ and $$\ell=1,2,3$$ (from left to right) as they are used in the localized corrector problem (5.1). Theorem 4.1 Let $$T\in {\mathcal{T}}_H$$ and $$v_H\in V_H$$ and let $${\mathcal{C}}_{T} v_H$$ denote the corresponding local subscale corrector as defined in (4.2). Then, we have |CTvH|H1(Ω∖ST,ℓ,b)≲βℓ|CTvH|H1(Ω). (4.5) The constant $$\beta$$ reads β=(4CIH(Hh)+3CIH(Hh)21+4CIH(Hh)+3CIH(Hh)2)1/2<1 (4.6) and is bounded away from 1. Before going to the proof of this theorem, we express the exponential decay in terms of patches in the following corollary. This is a direct consequence of Theorem 4.1 and the definition of $${\it{\Omega}}_{T,\ell,b}$$. Corollary 4.2 Let $$T\in {\mathcal{T}}_H$$ and $$v_H\in V_H$$ and let $${\mathcal{C}}_{T} v_H$$ denote the corresponding local subscale corrector as defined in (4.2). Then, we have |CTvH|H1(Ω∖ΩT,ℓ,b)≤|CTvH|H1(Ω∖ST,ℓ,b)≲βℓ|CTvH|H1(Ω) (4.7) with $$\beta<1$$ from (4.6). Remark 4.3 Recall that in this two-dimensional situation, $${C_{I_H}\big(\frac{H}{h}\big)}\lesssim \log(\tfrac{H}{h})$$ for $${C_{I_H}\big(\frac{H}{h}\big)}$$ from (3.1). Then given fixed $$H\in (0,1)$$ and let $$h\to 0$$, the constant $$\beta$$ scales as 1−β2≳11+log(h)2. In the three-dimensional case, Theorem 4.1 could essentially be proven in the same way, but the dependence of $${C_{I_H}\big(\frac{H}{h}\big)}$$ on $$H/h$$ is algebraic, so that the decay rate deteriorates very fast. Proof of Theorem 4.1. The crucial point in the proof is (4.10) below, which exploits the direction of $$b$$. This allows for patches that are only enlarged in the direction of $$-b$$. The remaining part of the proof then essentially follows as in Målqvist & Peterseim (2011). Define a cutoff function η:=1−η1η2, where $$0\leq\eta_1(x)\leq 1$$ and $$0\leq\eta_2(x)\leq 1$$ are one-dimensional continuous piecewise affine cutoff functions along $$t$$ and $$b$$, respectively. Recall that $$\mathrm{mid}_T$$ denotes the barycenter of a coarse element $$T$$, $$|b|=1$$ and $$t$$ is a unit vector orthogonal to $$b$$. We define $$\eta_1$$ and $$\eta_2$$ by η1(x)={1 if |(x−midT)⋅t|≤(ℓ−1)H,0 if |(x−midT)⋅t|≥ℓH (4.8) and η2(x)={1 if −(ℓ−1)H≤−(x−midT)⋅b≤(ℓ−1)H2ϵ,0 if −(x−midT)⋅b≥ℓH2ϵ or −(x−midT)⋅b≤−ℓH. (4.9) We obtain from the construction above that $$\nabla\eta_1(x)\cdot b=0$$ for all $$x\in{\it{\Omega}}$$ and $$\eta_1\leq 1$$. Moreover, because $$-(b\cdot \nabla \eta_2(x))\leq 0$$ if $$0\leq (x-\mathrm{mid}_T)\cdot b$$, we deduce −b⋅∇η=−(b⋅∇η1)η2−(b⋅∇η2)η1=−(b⋅∇η2)η1≤ϵH2. (4.10) Furthermore, $$\eta\vert_{S_{T,\ell-1,b}}=0$$ and $$\eta\vert_{{\it{\Omega}}\setminus S_{T,\ell,b}}=1$$, and $$\eta$$ is bounded between 0 and 1 and satisfies the Lipschitz continuity ‖∇η‖L∞(Ω)≤2/H. (4.11) Note that $$\mathrm{supp}(\nabla\eta)\subset S_{T,\ell,b}\setminus S_{T,\ell-1,b}$$. Let $$(\bullet,\bullet):=(\bullet,\bullet)_{L^2({\it{\Omega}})}$$ denote the $$L^2$$ scalar product and define $$\varphi:={\mathcal{C}}_{T} v_H$$. As $$\nabla\cdot b=0$$, we have $$(b\cdot\nabla(\eta\varphi),\eta\varphi)=0$$, and ϵ|φ|H1(Ω∖ST,ℓ,b)2 ≤ϵ(∇(ηφ),∇(ηφ))+(b⋅∇(ηφ),ηφ) =ϵ(∇φ,η∇(ηφ))+ϵ(∇η,φ∇(ηφ))+(b⋅∇(ηφ),ηφ) =ϵ(∇φ,∇(η2φ))+(b⋅∇(η2φ),φ)−ϵ(∇φ,ηφ∇η) +ϵ(∇η,φ∇(ηφ))−(b⋅∇η,ηφ2). Observe that $$\eta^2\varphi\in \mathrm{Ker}\,I_H$$, and we obtain ϵ(∇φ,∇(η2φ))+(b⋅∇(η2φ),φ)=a(η2φ,φ)=aT(η2φ,vH)=0 (4.12) by the definition of $${\mathcal{C}}_{T} $$ in (4.2). Thus, we arrive at ϵ|φ|H1(Ω∖ST,ℓ,b)2≤ϵ|(∇φ,ηφ∇η)|+ϵ|(∇η,φ∇(ηφ))|−(b⋅∇η,ηφ2). (4.13) We will estimate each term on the right-hand side of (4.13). With $$\eta\leq 1$$ and (4.11), a Cauchy inequality leads to ϵ|(∇φ,ηφ∇η)| ≤2ϵH−1|φ|H1(ST,ℓ,b∖ST,ℓ−1,b)‖φ‖L2(ST,ℓ,b∖ST,ℓ−1,b) ≤2CIH(Hh)ϵ|φ|H1(ST,ℓ,b∖ST,ℓ−1,b)2, where we have used the fact that $$\varphi\in \mathrm{Ker}\,I_H$$ and estimate (3.1) in the last inequality. The same arguments imply for the second term in (4.13), ϵ|(∇η,φ∇(ηφ))| ≤2ϵH−1‖φ‖L2(ST,ℓ,b∖ST,ℓ−1,b)|ηφ|H1(ST,ℓ,b∖ST,ℓ−1,b) ≤2ϵ(CIH(Hh)2+CIH(Hh))|φ|H1(ST,ℓ,b∖ST,ℓ−1,b)2. The crucial point in the estimation of the last term in (4.13) is estimate (4.10), which implies together with $$\eta\varphi^2\geq 0$$, −(b⋅∇η,ηφ2) ≤ϵH2‖φ‖L2(ST,ℓ,b∖ST,ℓ−1,b)2 ≤ϵCIH(Hh)2|φ|H1(ST,ℓ,b∖ST,ℓ−1,b)2. Assemble all estimates above for (4.13) to conclude ϵ|φ|H1(Ω∖ST,ℓ,b)2≤ϵ(4CIH(Hh)+3CIH(Hh)2)|φ|H1(ST,ℓ,b∖ST,ℓ−1,b)2. Define $$C\hspace{-0.6ex}\left(\tfrac{H}{h}\right):=4{C_{I_H}\big(\frac{H}{h}\big)}+3{C_{I_H}\big(\frac{H}{h}\big)}^2$$, which leads to |φ|H1(Ω∖ST,ℓ,b)2≤C(Hh)(|φ|H1(Ω∖ST,ℓ−1,b)2−|φ|H1(Ω∖ST,ℓ,b)2), and therefore |φ|H1(Ω∖ST,ℓ,b)2≤C(Hh)1+C(Hh)|φ|H1(Ω∖ST,ℓ−1,b)2. Repeating this process, we arrive at |φ|H1(Ω∖ST,ℓ,b)2≤(C(Hh)1+C(Hh))ℓ|φ|H1(Ω)2. This concludes the proof. □ Remark 4.4 (Nonconstant $$b$$) If the velocity field $$b$$ is divergence-free, but not globally constant, the definition of the rectangles $$S_{T,\ell,b}$$ has to be modified in that they have to follow the velocity. To avoid the computation of those patches for a complicated $$b$$, one may alternatively enlarge the patches adaptively using a posteriori error estimators (Larson & Målqvist, 2009). The following sketch shows that the proof of the exponential decay can be generalized to the situation where there exists a bounded diffeomorphism with bounded inverse that maps a constant reference velocity field $$b_\mathrm{ref}$$ to $$b$$, in the following sense. Assume that there exists a reference domain $${\it{\Omega}}_\mathrm{ref}$$ and a diffeomorphism $$\psi:{\it{\Omega}}_{\mathrm{ref}}\to {\it{\Omega}}$$, $$\psi\in \mathcal{C}^1({\it{\Omega}}_{\mathrm{ref}})$$, such that Dψ(y)bref=b(ψ(y))for all y∈Ωref, i.e., $$\psi$$ follows $$b$$ in the following sense. For any $$y_{0,\mathrm{ref}}\in{\it{\Omega}}_\mathrm{ref}$$, the function $$y(t):=\psi(y_{0,\mathrm{ref}}+t b_\mathrm{ref})$$ solves the ordinary differential equation (ODE) $$y'(t)=b(y(t))$$, i.e., $$y$$ follows $$b$$. The domain $$S_{T,\ell,b}$$ (formerly a rectangle) is then defined as $$S_{T,\ell,b}:=\psi(S_{\mathrm{ref},T,\ell,b_\mathrm{ref}})$$, where $$S_{\mathrm{ref},T,\ell,b_\mathrm{ref}}\subset {\it{\Omega}}_\mathrm{ref}$$ is defined for the constant vector field $$b_\mathrm{ref}$$ as in (4.3). The cutoff function $$\eta=1-\eta_1 \eta_2$$ is then defined by $$\eta_j(x):=\eta_{j,\mathrm{ref}}(\psi^{-1}(x))$$ for $$\eta_{j,\mathrm{ref}}$$ defined as in (4.8)–(4.9). The boundedness of $$D\psi^{-1}$$ then proves ‖∇η‖L∞(Ω)≲H−1. The definitions of $$\eta_1$$ and $$\eta_2$$ lead for all $$x\in{\it{\Omega}}$$ to b(x)⋅∇ηj(x)=∇ηj,ref|ψ−1(x)⋅(Dψ−1(x)b(x))=∇ηj,ref|ψ−1(x)⋅bref, which implies $$-b\cdot\nabla\eta\leq \epsilon/H^2$$. Theorem 4.1 then follows as before. Figure 4 displays the modified basis functions $$(1-{\mathcal{C}})\lambda_z$$ for $$z=(0.875,0.5)$$ for the following nonconstant vector fields b1(x)=5(x2−0.50.5−x1)andb2(x)=(2(⌈x−(0.50.5)⌉−1)mod2)b1(x), (4.14) where $$\lceil r\rceil:=\min\{k\in\mathbb{N}\mid k\geq r\}$$ denotes the ceiling function. For the first example, there exists a diffeomorphism $$\psi$$ as above in a subdomain $$(0.6,1)\times (0,1)$$, whereas for the second example, this is not true due to the jumps of $$b_2$$. Nevertheless, one observes a typical exponential decay pattern that follows $$b$$ very closely. Fig. 4. View largeDownload slide Top view of the modulus of (ideal) test basis functions $$(1-{\mathcal{C}})\lambda_z$$ with a logarithmic color scale. The underlying data are $$\epsilon=2^{-8}$$ and $$b_1$$ (left) and $$b_2$$ (right) from (4.14). Fig. 4. View largeDownload slide Top view of the modulus of (ideal) test basis functions $$(1-{\mathcal{C}})\lambda_z$$ with a logarithmic color scale. The underlying data are $$\epsilon=2^{-8}$$ and $$b_1$$ (left) and $$b_2$$ (right) from (4.14). 5. LOD method and error analysis On the basis of the results above, we conclude that the energy norm of $${\mathcal{C}}_T v$$ decreases very fast outside a local region around $$T$$ for any $$v\in V_H$$. Therefore, a localization process is feasible to reduce the computational costs of the ideal method but maintain a good accuracy. In this section, we want to localize the corrector problems (3.2). To this end, instead of solving them on the global domain $${\it{\Omega}}$$, we obtain a good approximation of those correctors by solving a local problem on $${\it{\Omega}}_{T,\ell,b}$$. First, let us introduce some notation. In the following, we will denote $$R_H=\mathrm{Ker}\,I_H$$ and $$R_H({\it{\Omega}}_{T,\ell,b})=\{w\in R_H,\text{ and }w=0 \text{ in }{\it{\Omega}}\setminus {\it{\Omega}}_{T,\ell,b}\}$$. Recall the local bilinear form $$a_\omega$$ defined in (4.1). The localized element corrector $${\mathcal{C}}_{T,\ell} : V_H\rightarrow R_H({\it{\Omega}}_{T,\ell,b})$$ is defined as follows: given $$v_H\in V_H$$, let $${\mathcal{C}}_{T,\ell}v_H\in R_H({\it{\Omega}}_{T,\ell,b})$$ satisfy aΩT,ℓ,b(w,CT,ℓvH)=aT(w,vH)for all w∈RH(ΩT,ℓ,b). (5.1) Then we denote $${\mathcal{C}}_{\ell} :=\sum\nolimits_{T\in {\mathcal{T}}_H}{\mathcal{C}}_{T,\ell}$$; see Fig. 5 for an illustration of the localized correctors $${\mathcal{C}}_{\ell} \lambda_z$$ and the corresponding localized test basis. Fig. 5. View largeDownload slide Localized element correctors $${\mathcal{C}}_{T,\ell}\lambda_z$$ for $$\ell=2$$ and all four elements $$T$$ adjacent to the vertex $$z=[0.5,0.5]$$ (top), localized nodal corrector $${\mathcal{C}}_{\ell}\lambda_z=\sum_{T\ni z} {\mathcal{C}}_{T,\ell}\lambda_z$$ (bottom left) and corresponding test basis function $$(1-{\mathcal{C}}_{\ell})\lambda_z$$ (bottom right). The underlying data are $$b=[\cos(0.7),\sin(0.7)]$$ and $$\epsilon=2^{-7}$$. The computations have been performed by standard linear finite elements on local fine meshes of width $$h=2^{-8}$$. See Fig. 1 for a comparison with the ideal global corrector and basis. Fig. 5. View largeDownload slide Localized element correctors $${\mathcal{C}}_{T,\ell}\lambda_z$$ for $$\ell=2$$ and all four elements $$T$$ adjacent to the vertex $$z=[0.5,0.5]$$ (top), localized nodal corrector $${\mathcal{C}}_{\ell}\lambda_z=\sum_{T\ni z} {\mathcal{C}}_{T,\ell}\lambda_z$$ (bottom left) and corresponding test basis function $$(1-{\mathcal{C}}_{\ell})\lambda_z$$ (bottom right). The underlying data are $$b=[\cos(0.7),\sin(0.7)]$$ and $$\epsilon=2^{-7}$$. The computations have been performed by standard linear finite elements on local fine meshes of width $$h=2^{-8}$$. See Fig. 1 for a comparison with the ideal global corrector and basis. In the following lemma, we will show that $${\mathcal{C}}_{T,\ell}$$ is a good approximation of $${\mathcal{C}}_{T} $$ provided that the local patches $${\it{\Omega}}_{T,\ell,b}$$ are sufficiently large. For ease of presentation, we denote the mesh Péclet number $${\mathrm{Pe}_{H,b,\epsilon}\,}$$ of $${\mathcal{T}}_H$$ by PeH,b,ϵ:=H‖b‖L∞(Ω)/ϵ. (5.2) Recall the definition of $$\beta$$ from (4.6). Lemma 5.1 Given $$v\in V_H$$ and $$\ell\in\mathbb{N}_{+}$$, it holds that |CTv−CT,ℓv|H1(Ω)≲(1+PeH,b,ϵCIH(Hh))2(CIH(Hh)+1)βℓ−1|v|H1(T). (5.3) Proof. Define $$e_{T,\ell}:={\mathcal{C}}_{T} v-{\mathcal{C}}_{T,\ell} v$$. In view of $$R_{H}({\it{\Omega}}_{T,\ell,b})\subset R_{H}$$, the definitions of the correctors in (5.1) and (3.2) and the orthogonality of Petrov–Galerkin type lead to ϵ|eT,ℓ|H1(Ω)2=a(eT,ℓ,eT,ℓ)=a(eT,ℓ−w,eT,ℓ) for all w∈RH(ΩT,ℓ,b). As $$I_H(e_{T,\ell})=0$$, Hölder’s inequality and the approximation property (3.1) of $$I_H$$ imply |eT,ℓ|H1(Ω)2≤(1+PeH,b,ϵCIH(Hh))|eT,ℓ|H1(Ω)|eT,ℓ−w|H1(Ω). As $$w\in R_{H}({\it{\Omega}}_{T,\ell,b})$$ is arbitrary, we arrive at |eT,ℓ|H1(Ω)≤(1+PeH,b,ϵCIH(Hh))|CTv−w|H1(Ω). (5.4) In the following, we construct a specific $$w\in R_H({\it{\Omega}}_{T,\ell,b})$$ to control the term $$\left|{{\mathcal{C}}_{T}v-w}\right|_{H^{1}\left({{\it{\Omega}}}\right)}$$. Let $$\eta$$ denote the cutoff function from the proof of Theorem 4.1, such that $$\eta\vert_{S_{T,\ell-1,b}}=0$$ and $$\eta\vert_{{\it{\Omega}}\setminus S_{T,\ell,b}}=1$$. Note that $$S_{T,\ell,b}\subset{\it{\Omega}}_{T,\ell,b}$$ and therefore $$\mu:=1-\eta$$ satisfies $$\mu\vert_{{\it{\Omega}}\setminus{\it{\Omega}}_{T,\ell,b}}=0$$. In addition, $$\mu$$ is bounded between 0 and 1 and satisfies the Lipschitz continuity ‖∇μ‖L∞(Ω)≤2H−1. (5.5) Define $$w=\mu{\mathcal{C}}_{T}v$$; then $$w\in R_H({\it{\Omega}}_{T,\ell,b})$$. As $${\mathcal{C}}_{T}v\in R_H$$, the fact that $$0\leq\mu\leq 1$$ and (5.5) lead as in the proof of Theorem 4.1 to |CTv−w|H1(Ω) =|CTv−μCTv|H1(Ω∖ST,ℓ−1,b) ≤2(CIH(Hh)+1)|CTv|H1(Ω∖ST,ℓ−1,b). Theorem 4.1 then implies |CTv−w|H1(Ω) ≲(CIH(Hh)+1)βℓ−1|CTv|H1(Ω). The combination with (5.4) implies |CTv−CT,ℓv|H1(Ω)≲(1+PeH,b,ϵCIH(Hh))(CIH(Hh)+1)βℓ−1|CTv|H1(Ω). In the end, we show the stability of $$\mathcal{C}_{T}$$ to bound the term $$\left|{{\mathcal{C}}_{T}v}\right|_{H^{1}\left({{\it{\Omega}}}\right)}$$. As $$I_H({\mathcal{C}}_{T}v)=0$$, the stability of $${\mathcal{C}}_{T}$$ follows from ϵ|CTv|H1(Ω)2 =a(CTv,CTv)=aT(CTv,v) ≤ϵ|CTv|H1(Ω)|v|H1(T)+‖b‖L∞(T)|v|H1(T)‖CTv‖L2(Ω) ≤(ϵ+H‖b‖L∞(T)CIH(Hh))|CTv|H1(Ω)|v|H1(T), where the definition of the element corrector in (5.1) implies the second equality and the approximation property (3.1) leads to the last inequality. This proves the assertion. □ The following theorem assembles the local estimates from Lemma 5.1 to derive an estimate for the global corrector. Theorem 5.2 Given $$v\in V_H$$ and $$\ell\in\mathbb{N}_{+}$$, it holds that |Cv−Cℓv|H1(Ω)≲C(H,h,ϵ,b,ℓ)βℓ−1|v|H1(Ω) (5.6) with C(H,h,ϵ,b,ℓ) :=(1+PeH,b,ϵCIH(Hh))2(CIH(Hh)+1) ×(1+2CIH(Hh)+PeH,b,ϵCIH(Hh))Col,ℓ+2(ϵ)1/2. (5.7) Proof. Set $$z:={\mathcal{C}} v-{\mathcal{C}}_{\ell-2}v\in \mathrm{Ker}\,I_H$$ and $$z_T:={\mathcal{C}}_{T}v-{\mathcal{C}}_{T,\ell-2}v$$, then $$z=\sum\limits_{T\in{\mathcal{T}}_H}z_T$$. We have ϵ|z|H1(Ω)2=∑T∈THa(z,zT). (5.8) We estimate $$a(z,z_T)$$ for each coarse element $$T\in \mathcal{T}_H$$. Recall that we defined a cutoff function $$\eta$$ in the proof of Theorem 4.1. Note that $${\it{\Omega}}_{T,\ell-2,b}\subset S_{T,\ell-1,b}$$. By construction, we have $$\eta z\in R_H({\it{\Omega}}\setminus S_{T,\ell-1,b})\subset R_H({\it{\Omega}}\setminus {\it{\Omega}}_{T,\ell-2,b})$$. As $${\mathcal{C}}_{T,\ell-2} v\vert_{{\it{\Omega}}\setminus{\it{\Omega}}_{T,\ell-2,b}}=0$$, this implies a(ηz,zT)=a(ηz,CTv). Furthermore, notice that $$\eta z\in \mathrm{ker}(I_H)$$, which combined with (4.2) yields a(ηz,CTv)=aT(ηz,v)=0. As a consequence, we obtain a(z,zT)=a(ηz,zT)+a((1−η)z,zT)=a((1−η)z,zT). In the following, we will bound the term $$a((1-\eta)z,z_T)$$. Recall from the proof of Theorem 4.1 that $$(1-\eta)\vert_{{\it{\Omega}}\setminus S_{T,\ell,b}}=0$$, $$\|\nabla(1-\eta)\|_{L^\infty({\it{\Omega}})}\leq 2H^{-1}$$ and $$\|(1-\eta)\|_{L^\infty({\it{\Omega}})}\leq 1$$. Taking into account that $$I_H(z)=I_H(z_T)=0$$, the stability of the projector $$I_H$$ from (3.1), therefore, leads to a((1−η)z,zT)≤ϵ|(1−η)z|H1(ST,ℓ,b)|zT|H1(ST,ℓ,b)+‖b‖L∞(Ω)‖z‖L2(ST,ℓ,b)|zT|H1(ST,ℓ,b) ≤(ϵ(1+2CIH(Hh))+‖b‖L∞(Ω)HCIH(Hh))|z|H1(ST,ℓ,b)|zT|H1(ST,ℓ,b). As $$S_{T,\ell,b}\subset {\it{\Omega}}_{T,\ell,b}$$, the combination with (5.8) and the application of a discrete Cauchy–Schwarz inequality yields |z|H1(Ω)2 ≤(1+2CIH(Hh)+PeH,b,ϵCIH(Hh))∑T∈TH|z|H1(ST,ℓ,b)|zT|H1(ST,ℓ,b) ≤(1+2CIH(Hh)+PeH,b,ϵCIH(Hh)) ×(∑T∈TH|z|H1(ΩT,ℓ,b)2)1/2(∑T∈TH|zT|H1(ST,ℓ,b)2)1/2. Lemma 5.1 implies (∑T∈TH|zT|H1(ST,ℓ,b)2)1/2≲(1+PeH,b,ϵCIH(Hh))2(CIH(Hh)+1)βℓ−3|v|H1(Ω), while the bounded overlap of the patches from (4.4) implies (∑T∈TH|z|H1(ΩT,ℓ,b)2)1/2≤Col,ℓ(ϵ)1/2|z|H1(Ω). In the end, the combination of the previously displayed inequalities and the shift $$\ell\mapsto \ell+2$$ shows the assertion. □ Now we are ready to define the localized multiscale test space as WH,ℓ=(1−Cℓ)VH. The Petrov–Galerkin method for the approximation of (2.5) based on the trial–test pairing $$(V_{H},W_{H,\ell})$$ defined above seeks $$u_{H,\ell}\in V_{H}$$ satisfying a(uH,ℓ,wH,ℓ)=⟨f,wH,ℓ⟩H−1(Ω)×H01(Ω)for all wH,ℓ∈WH,ℓ. (5.9) Lemma 5.3 (Inf–sup stability) If $$\ell$$ is sufficiently large, i.e., the oversampling condition ℓ≳1+|log(CIH(Hh))|+|log(C(H,h,ϵ,b,ℓ))|+|log(1+CIH(Hh)+PeH,b,ϵCIH(Hh)2)||log(4CIH(Hh)+3CIH(Hh)2)−log(1+4CIH(Hh)+3CIH(Hh)2)| (5.10) is satisfied, then the Petrov–Galerkin method (5.9) is inf–sup stable and infwH,ℓ∈WH,ℓ∖{0}supuH∈VH∖{0}a(uH,wH,ℓ)|uH|H1(Ω)|wH,ℓ|H1(Ω)≳ϵCIH(Hh). (5.11) Remark 5.4 If $$H/h=1/\epsilon$$ and $$\lvert b\rvert=1$$ then (5.10) reads ℓ≳(log(ϵ))2, i.e., the local patch size $$\ell$$ depends on $$\log(\epsilon)$$ algebraically. Remark 5.5 As the dimension of $$V_H$$ equals the dimension of $$W_{H,\ell}$$, the reverse inf–sup condition infuH∈VH∖{0}supwH,ℓ∈WH,ℓ∖{0}a(uH,wH,ℓ)|uH|H1(Ω)|wH,ℓ|H1(Ω)≳ϵCIH(Hh) (5.12) follows from Lemma 5.3. Proof of Lemma 5.3. Let $$w_{H,\ell}\in W_{H,\ell}$$, and set $$w_{H}=(1-{\mathcal{C}} )I_{H}w_{H,\ell}\in W_H$$. By Lemma 3.4, there exists $$u_H\in V_H$$, such that a(uH,wH)≥ϵCIH(Hh)|uH|H1(Ω)|wH|H1(Ω). (5.13) Taking into account that $$w_{H,\ell}=I_H w_{H,\ell}-{\mathcal{C}}_\ell I_H w_{H,\ell}$$, we arrive at $$w_H-w_{H,\ell} = ({\mathcal{C}}_\ell-{\mathcal{C}})I_H w_{H,\ell}$$. As a consequence, Theorem 5.2 together with the stability of $$I_H$$ from (3.1) implies |wH−wH,ℓ|H1(Ω) ≤C~C(H,h,ϵ,b,ℓ)βℓ−1|IHwH,ℓ|H1(Ω) ≤C~C(H,h,ϵ,b,ℓ)CIH(Hh)βℓ−1|wH,ℓ|H1(Ω). Here, $$\tilde{C}$$ denotes the constant hidden in $$\lesssim$$ in Theorem 5.2, which is independent of $$H$$, $$h$$ or $$\epsilon$$. The combination with a triangle inequality leads to |wH|H1(Ω) ≥|wH,ℓ|H1(Ω)−|wH−wH,ℓ|H1(Ω) ≥(1−C~C(H,h,ϵ,b,ℓ)CIH(Hh)βℓ−1)|wH,ℓ|H1(Ω). As $$I_H(w_{H,\ell}-w_H)=0$$, i.e., $$w_{H,\ell}-w_H\in R_H$$, this leads to |a(uH,wH,ℓ−wH)| ≤(ϵ+‖b‖L∞(Ω)HCIH(Hh))|uH|H1(Ω)|wH,ℓ−wH|H1(Ω) =ϵ(1+PeH,b,ϵCIH(Hh))|uH|H1(Ω)|wH,ℓ−wH|H1(Ω). The combination of the above displayed inequalities results in a(uH,wH,ℓ) =a(uH,wH)+a(uH,wH,ℓ−wH) ≥ϵCIH(Hh)(1−C~C(H,h,ϵ,b,ℓ)CIH(Hh)βℓ−1)|uH|H1(Ω)|wH,ℓ|H1(Ω) −ϵ(1+PeH,b,ϵCIH(Hh))C(H,h,ϵ,b,ℓ)C~CIH(Hh)βℓ−1|uH|H1(Ω)|wH,ℓ|H1(Ω). Recall the definition of $$\beta$$ from (4.6). If $$\ell$$ satisfies (5.10) then we obtain (5.11). □ We are ready to estimate the error $$\left|{u_{H}-u_{H,\ell}}\right|_{H^{1}\left({{\it{\Omega}}}\right)}$$ coming from the localization. Lemma 5.6 Let $$\ell$$ satisfy (5.10). Then |uH−uH,ℓ|H1(Ω) ≲CIH(Hh)2C(H,h,ϵ,b,ℓ)(1+PeH,b,ϵCIH(Hh)) ×βℓ−1|uh−uH|H1(Ω). (5.14) Proof. Notice that $$u_{H}-u_{H,\ell}\in V_H$$ is a coarse finite element function. Therefore, the inf–sup condition (5.12) guarantees the existence of $$w_{H,\ell}\in W_{H,\ell}$$ with |uH−uH,ℓ|H1(Ω)≲CIH(Hh)ϵa(uH−uH,ℓ,wH,ℓ)|wH,ℓ|H1(Ω). In view of $$w_{H,\ell}\in W_{H,\ell}\subset V_h$$, the standard Galerkin problem (2.5) and the VMS method (5.9) imply a(uH−uH,ℓ,wH,ℓ)=a(uH−uh,wH,ℓ). Define $$w_{H}:=I_H w_{H,\ell} - {\mathcal{C}}(I_{H}w_{H,\ell})\in W_{H}\subset V_h$$. Together with the orthogonality of Petrov–Galerkin type, we obtain a(uH−uh,wH,ℓ)=a(uH−uh,wH,ℓ−wH). Taking into account that $$w_{H,\ell}-w_H= ({\mathcal{C}}-{\mathcal{C}}_\ell) I_H w_{H,\ell}$$, the combination with a Cauchy inequality, $$w_{H,\ell}-w_H\in\mathrm{ker}(I_H)$$ and an application of Theorem 5.2 lead to a(uH−uh,wH,ℓ−wH)≲(ϵ+‖b‖L∞(Ω)HCIH(Hh))C(H,h,ϵ,b,ℓ)βℓ−1|IHwH,ℓ|H1(Ω)|uh−uH|H1(Ω). The stability of $$I_H$$ from (3.1) implies the assertion. □ Lemma 5.6 allows us to bound the error for the localized VMS method in the following manner. Theorem 5.7 (Global error estimate for the localized VMS method) Let $$\ell$$ satisfy (5.10); then |uh−uH,ℓ|H1(Ω) ≲(CIH(Hh)+CIH(Hh)3C(H,h,ϵ,b,ℓ)(1+PeH,b,ϵCIH(Hh))βℓ−1) ×minvH∈VH|uh−vH|H1(Ω) with the constant $$C(H,h,\epsilon,b,\ell)$$ from (5.7). Proof. The proof follows directly from a triangle inequality, Proposition 3.1 and Lemma 5.6. □ Although Theorem 5.7 provides a best approximation result, the assertion still depends on $$\epsilon$$, which is hidden in the best approximation $$\min_{v_H\in V_H} \left|{u_h-v_H}\right|_{H^{1}\left({{\it{\Omega}}}\right)}$$. The locality in the error bound of the ideal method from Proposition 3.1 transfers to the VMS method defined in (5.9) and results in the local error bound in the following theorem. Note that the error from the localization still depends on the mesh Péclet number of $${\mathcal{T}}_H$$ and still contains the best approximation error on the whole domain. Nevertheless, these ill-behaved terms are weighted by the exponentially decaying term $$\beta^{\ell-1}$$, where $$\beta$$ is bounded above from 1. Theorem 5.8 (Local error estimate for the localized VMS method) Let $$\ell$$ satisfy (5.10). Then for any $$\mathcal{K}\subset {\mathcal{T}}_H$$ and $$\omega:=\cup\mathcal{K}$$, it holds that |uh−uH,ℓ|H1(ω) ≲CIH(Hh)minvH∈VH|uh−vH|H1(ω) +CIH(Hh)3C(H,h,ϵ,b,ℓ)(1+PeH,b,ϵCIH(Hh))βℓ−1minvH∈VH|uh−vH|H1(Ω). Remark 5.9 (Complexity) Problem (5.9) on the coarse scale consists of $$\mathcal{O}(1/H^2)$$ degrees of freedom (DOFs). Corresponding to each of those DOFs, one localized corrector problem (5.1) has to be solved, which relates to $$\mathcal{O}(\ell^2 H^3/ (h^2\epsilon))$$ DOFs in the worst-case scenario. If the mesh is structured, the number of corrector problems that have to be solved can be reduced to $$\mathcal{O}(\ell H/\epsilon)$$ (cf. Gallistl & Peterseim, 2015). 6. Numerical experiment In this section, we present one simple numerical test to illustrate the theoretical convergence results of the localized method proposed in (5.9). We take $${\it{\Omega}}=(0,1)\times(0,1)$$, the velocity field $$b=(\cos(0.7), \sin(0.7))^{\rm T}$$, the volume force $$f\equiv 1$$ and $$\epsilon = 2^{-7}$$. The reference solution $$u_h$$ is obtained through (2.5) by taking $$h=\sqrt{2}\,2^{-8}$$. We will compare our approach with the SUPG method. Let us briefly review the SUPG model to (1.1) (Franca et al., 1992). Let $$(\bullet,\bullet)_T:=(\bullet,\bullet)_{L^2(T)}$$ denote the $$L^2$$ scalar product over a triangle $$T\in{\mathcal{T}}_H$$. Then the SUPG method seeks $$u_H\in V_H$$ such that BSUPG(uHSUPG,vH)=FSUPG(vH)for all vH∈VH (6.1) with BSUPG(uHSUPG,vH)=a(uHSUPG,vH)+δSUPG∑T∈TH(b⋅∇uHSUPG,b⋅∇vH)T and FSUPG(vH)=⟨f,vH⟩H−1(Ω)×H01(Ω)+δSUPG∑k∈TH(f,b⋅∇vH)T. Here, $$\delta_\mathrm{SUPG}$$ indicates the stability parameter, and we choose δSUPG=H8max(ϵ,H/2) in our numerical test. The reference solution from (2.5) and the coarse-scale solution from (5.9) and the SUPG solution from (6.1) with $$H=\sqrt{2}\,2^{-4}$$ are depicted in Fig. 6. One can observe that the classical FEM approximation with $$H=\sqrt{2}\,2^{-4}$$ is not stable around the boundary layers (i.e., the top and right boundaries) and shows spurious oscillations and thus fails to provide a reliable solution. Nevertheless, both the SUPG method and the ideal method are stable and generate an accurate solution. In Fig. 7, we display the solutions for fixed $$y=0.75$$ to illustrate the stability and accuracy of the VMS method. We observe large oscillations in the coarse-scale solution obtained through classical FEM when $$x$$ approaches 1, while the SUPG and the VMS methods yield reliable solutions. The smearing is restricted to one layer of elements around the boundary. We can also conclude that the SUPG and the VMS methods reproduce the reference solution away from $$x=1$$ and the latter shows slightly less smearing. We want to highlight that the localization parameter is $$\ell=1$$ for the VMS method in this example. Fig. 6. View largeDownload slide Reference solution (top left), classical FEM approximation (top right), SUPG approximation (bottom left) and multiscale approximation for $$\ell=1$$ (bottom right) for $$\epsilon=2^{-7}$$ and $$H=\sqrt{2}\,2^{-4}$$. Fig. 6. View largeDownload slide Reference solution (top left), classical FEM approximation (top right), SUPG approximation (bottom left) and multiscale approximation for $$\ell=1$$ (bottom right) for $$\epsilon=2^{-7}$$ and $$H=\sqrt{2}\,2^{-4}$$. Fig. 7. View largeDownload slide Reference solution, classical FEM, SUPG and multiscale approximation for $$\ell=1$$ at the line $$y=0.75$$ on a mesh with mesh size $$H=\sqrt{2}\,2^{-4}$$ for $$\epsilon=2^{-7}$$. Fig. 7. View largeDownload slide Reference solution, classical FEM, SUPG and multiscale approximation for $$\ell=1$$ at the line $$y=0.75$$ on a mesh with mesh size $$H=\sqrt{2}\,2^{-4}$$ for $$\epsilon=2^{-7}$$. Tables 1 and 2 list the errors between the localized solutions (5.9) and the reference solution $$u_h$$ under various coarse mesh sizes $$H$$ and localization parameters $$\ell$$. We observe an optimal convergence rate of $$\mathcal{O}(H)$$ in Table 1 for the error in the $$H^1$$ seminorm in the domain $$[0,0.75]\times [0,0.75]$$ away from the boundary layers and an optimal convergence rate of $$\mathcal{O}(H^2)$$ in Table 2 for the global error in the $$L^2$$ norm. Although Theorem 5.8 guarantees optimality only under the assumption that $$\ell$$ is large enough in the sense of (5.10), the numerical experiment demonstrates that $$\ell=1$$ is sufficient for an accurate solution, which implies a potentially huge computational reduction. Table 1 The error $$\|\nabla(u_h-u_{H,\ell})\|_{L^2({\it{\Omega}}_r)}$$ for $${\it{\Omega}}_r=[0,0.75]\times [0,0.75]$$ for different localization parameters $$\ell$$ and mesh sizes $$H$$ for $$\epsilon=2^{-7}$$ $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $$\ell=5$$ $$\ell=6$$ $$H=\sqrt{2}\,2^{-3}$$ 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 $$H=\sqrt{2}\,2^{-4}$$ 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 $$H=\sqrt{2}\,2^{-5}$$ 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 $$H=\sqrt{2}\,2^{-6}$$ 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $$\ell=5$$ $$\ell=6$$ $$H=\sqrt{2}\,2^{-3}$$ 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 $$H=\sqrt{2}\,2^{-4}$$ 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 $$H=\sqrt{2}\,2^{-5}$$ 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 $$H=\sqrt{2}\,2^{-6}$$ 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 Table 1 The error $$\|\nabla(u_h-u_{H,\ell})\|_{L^2({\it{\Omega}}_r)}$$ for $${\it{\Omega}}_r=[0,0.75]\times [0,0.75]$$ for different localization parameters $$\ell$$ and mesh sizes $$H$$ for $$\epsilon=2^{-7}$$ $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $$\ell=5$$ $$\ell=6$$ $$H=\sqrt{2}\,2^{-3}$$ 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 $$H=\sqrt{2}\,2^{-4}$$ 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 $$H=\sqrt{2}\,2^{-5}$$ 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 $$H=\sqrt{2}\,2^{-6}$$ 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $$\ell=5$$ $$\ell=6$$ $$H=\sqrt{2}\,2^{-3}$$ 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 5.14e$$-$$02 $$H=\sqrt{2}\,2^{-4}$$ 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 2.57e$$-$$02 $$H=\sqrt{2}\,2^{-5}$$ 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 1.27e$$-$$02 $$H=\sqrt{2}\,2^{-6}$$ 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 6.23e$$-$$03 Table 2 The error $$\|u_h-u_{H,\ell}\|_{L^2({\it{\Omega}})}$$ for different localization parameters $$\ell$$ and mesh sizes $$H$$ for $$\epsilon=2^{-7}$$ $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $$\ell=5$$ $$\ell=6$$ $$H=0.17678$$ 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 $$H=0.088388$$ 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 $$H=0.044194$$ 2.31e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 $$H=0.022097$$ 7.25e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $$\ell=5$$ $$\ell=6$$ $$H=0.17678$$ 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 $$H=0.088388$$ 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 $$H=0.044194$$ 2.31e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 $$H=0.022097$$ 7.25e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 Table 2 The error $$\|u_h-u_{H,\ell}\|_{L^2({\it{\Omega}})}$$ for different localization parameters $$\ell$$ and mesh sizes $$H$$ for $$\epsilon=2^{-7}$$ $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $$\ell=5$$ $$\ell=6$$ $$H=0.17678$$ 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 $$H=0.088388$$ 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 $$H=0.044194$$ 2.31e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 $$H=0.022097$$ 7.25e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $$\ell=5$$ $$\ell=6$$ $$H=0.17678$$ 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 9.45e$$-$$02 $$H=0.088388$$ 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 5.34e$$-$$02 $$H=0.044194$$ 2.31e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 2.32e$$-$$02 $$H=0.022097$$ 7.25e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 7.27e$$-$$03 The convergence rate for $$u_{H,1}$$ with various $$\epsilon$$ in a range from $$2^{-5}$$ to $$2^{-8}$$ is shown in Figs 8 and 9. The error is stable and of order $$\mathcal{O}(H)$$ with respect to the $$H^1$$ seminorm in a region away from boundary layers and of order $$\mathcal{O}(H^2)$$ in the global $$L^2$$ norm with a preasymptotic effect for smaller values of $$\epsilon$$. For comparison, the nodal interpolation error (i.e., the error from the ideal method) in the global $$L^2$$ norm is depicted, which agrees with $$\left\|{u_h-u_{H,1}}\right\|_{L^2({\it{\Omega}})}$$ very well. This justifies the fast convergence of the localized method with respect to the localization parameter $$\ell$$ for all of the considered values of $$\epsilon$$. Fig. 8. View largeDownload slide The errors $$\|\nabla(u_h-u_{H,1})\|_{L^2({\it{\Omega}}_r)}$$ for $${\it{\Omega}}_r=[0,0.75]\times [0,0.75]$$ for different values of $$\epsilon$$. Fig. 8. View largeDownload slide The errors $$\|\nabla(u_h-u_{H,1})\|_{L^2({\it{\Omega}}_r)}$$ for $${\it{\Omega}}_r=[0,0.75]\times [0,0.75]$$ for different values of $$\epsilon$$. Fig. 9. View largeDownload slide The errors $$\|u_h-u_{H,1}\|_{L^2({\it{\Omega}})}$$ and $$\|u_h-I_H u_h\|_{L^2({\it{\Omega}})}$$ for different values of $$\epsilon$$. Fig. 9. View largeDownload slide The errors $$\|u_h-u_{H,1}\|_{L^2({\it{\Omega}})}$$ and $$\|u_h-I_H u_h\|_{L^2({\it{\Omega}})}$$ for different values of $$\epsilon$$. 7. Conclusions In this article, a singularly perturbed convection–diffusion equation was considered, and we obtained a stable locally quasi-optimal variational multiscale method based on the nodal interpolation operator. Because of the high complexity involved in solving the global correctors, which account for the main component of the VMS method, a further model reduction was proceeded by localization techniques based on the LOD method. This localization employs local patches that depend on the velocity field $$b$$ and the singular perturbation parameter $$\epsilon$$. The error of the localization decays exponentially. We also provided a numerical experiment to illustrate our theoretical results. The stability constant of the nodal interpolation operator that occurs in the error estimate depends logarithmically on $$H/h$$ (and so on $$\epsilon$$). In the three-dimensional case, this stability estimate depends polynomially on $$H/h$$. Therefore, a generalization of the proposed method to three dimensions does not seem reasonable. The local patches in the localized computation of the corrector depend on $$\epsilon$$. It is an open question whether this is optimal or whether a further reduction or simplification is possible. Acknowledgements This article has been written on the occasion of the trimester program on multiscale problems of the Hausdorff Institut for Mathematics (Bonn). Funding Hausdorff Center for Mathematics Bonn (to G.L.). Deutsche Forschungsgemeinschaft in the Priority Program 1748 ‘Reliable simulation techniques in solid mechanics: Development of nonstandard discretization methods, mechanical and mathematical analysis’ under the project ‘Adaptive isogeometric modeling of propagating strong discontinuities in heterogeneous materials’ (PE2143/2-1; to D.P.). References Barrenechea G. R. , Burman E. , & Karakatsani F. , ( 2016 ) Edge-based nonlinear diffusion for finite element approximations of convection–diffusion equations and its relation to algebraic flux-correction schemes. Numer. Math. , 135 , 521 – 545 . Google Scholar CrossRef Search ADS Brezzi F. , Marini D. , & Süli E. ( 2000 ) Residual-free bubbles for advection-diffusion problems: the general error analysis. Numer. Math. , 85 , 31 – 47 . Google Scholar CrossRef Search ADS Calo V. M. , Chung E. T. , Efendiev Y. , & Leung W. T. ( 2016 ) Multiscale stabilization for convection-dominated diffusion in heterogeneous media. Comput. Methods Appl. Mech. Eng. , 304 , 359 – 377 . Google Scholar CrossRef Search ADS Cangiani A. & Süli E. ( 2005a ) Enhanced RFB method. Numer. Math. , 101 , 273 – 308 . Google Scholar CrossRef Search ADS Cangiani A. & Süli E. ( 2005b ) Enhanced residual-free bubble method for convection-diffusion problems. Int. J. Numer. Methods Fluids , 47 , 1307 – 1313 . 8th ICFD Conference on Numerical Methods for Fluid Dynamics. Part 2 . Google Scholar CrossRef Search ADS Christiansen S. , Halvorsen T. , & Sørensen T. , ( 2016 ) Stability of an upwind Petrov-Galerkin discretization of convection diffusion equations. ArXiv e-prints, arXiv:1406.0390v2 . Codina R. ( 2000 ) Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods. Comput. Methods Appl. Mech. Eng. , 190 , 1579 – 1599 . Google Scholar CrossRef Search ADS Demkowicz L. , Gopalakrishnan J. , & Niemi A. H. ( 2012 ) A class of discontinuous Petrov-Galerkin methods. Part III: adaptivity. Appl. Numer. Math. , 62 , 396 – 427 . Google Scholar CrossRef Search ADS Elfverson D. ( 2015 ) A discontinuous Galerkin multiscale method for convection-diffustion problems. ArXiv e-prints . Franca L. P. , Frey S. L. & Hughes T. J. ( 1992 ) Stabilized finite element methods: I. Application to the advective-diffusive model. Comput. Methods Appl. Mech. Eng. , 95 , 253 – 276 . Google Scholar CrossRef Search ADS Gallistl D. & Peterseim D. , ( 2015 ) Stable multiscale Petrov-Galerkin finite element method for high frequency acoustic scattering. Comp. Meth. Appl. Mech. Eng. , 295 , 1 – 17 . Google Scholar CrossRef Search ADS Harder C. , Paredes D. , & Valentin F. , ( 2015 ) On a multiscale hybrid-mixed method for advective-reactive dominated problems with heterogeneous coefficients. Multiscale Model. Simul. , 13 , 491 – 518 . Google Scholar CrossRef Search ADS Hughes T. J. R. & Sangalli G. , ( 2007 ) Variational multiscale analysis: the fine-scale Green’s function, projection, optimization, localization, and stabilized methods. SIAM J. Numer. Anal. , 45 , 539 – 557 . Google Scholar CrossRef Search ADS John V. & Schmeyer E. , ( 2009 ) On finite element methods for 3D time dependent convection-diffusion-reaction equations with small diffusion. BAIL 2008 - boundary and interior layers , 173 – 181 , Lecture Notes in Computational Science and Engineering , 69 , Berlin : Springer . Google Scholar CrossRef Search ADS Larson M. G. & Måalqvist A. , ( 2009 ) An adaptive variational multiscale method for convection-diffusion problems. Comm. Numer. Methods Eng. , 25 , 65 – 79 . Google Scholar CrossRef Search ADS Måalqvist A. & Peterseim D. , ( 2014 ) Localization of elliptic multiscale problems. Math. Comp. , 83 , 2583 – 2603 . Google Scholar CrossRef Search ADS Melenk J. M. ( 1997 ) On the robust exponential convergence of hp finite element methods for problems with boundary layers. IMA J. Numer. Anal. , 17 , 577 – 601 . Google Scholar CrossRef Search ADS Melenk J. M. ( 2002 ) hp-Finite Element Methods for Singular Perturbations . Lecture Notes in Mathematics , vol. 1796 . Berlin : Springer . Google Scholar CrossRef Search ADS Park P. J. & Hou T. Y. ( 2004 ) Multiscale numerical methods for singularly perturbed convection-diffusion equations. Int. J. Comput. Methods , 1 , 17 – 65 . Google Scholar CrossRef Search ADS Peterseim D. ( 2017 ) Eliminating the pollution effect in Helmholtz problems by local subscale correction. Math. Comp. , 86 , 1005 – 1036 . Google Scholar CrossRef Search ADS Peterseim D. ( 2016b ) Variational multiscale stabilization and the exponential decay of fine-scale correctors. Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations ( Barrenechea G. R. , Brezzi F. , Cangiani A. , & Georgoulis E. H. eds). Lecture Notes in Computational Science and Engineering , vol. 114 . Cham: Springer , pp. 343 – 369 . Google Scholar CrossRef Search ADS Qiu W. & Shi K. , ( 2015 ) An HDG method for convection diffusion equation. J. Sci. Comput. , 66 , 346 – 357 . Google Scholar CrossRef Search ADS Roos H. , Stynes M. , & Tobiska L. , ( 2008 ) Robust Numerical Methods for Singularly Perturbed Differential Equations: Convection-Diffusion-Reaction and Flow Problems. Springer Series in Computational Mathematics , 24 . Berlin: Springer-Verlag . Yserentant H. ( 1986 ) On the multilevel splitting of finite element spaces. Numer. Math. , 49 , 379 – 412 . Google Scholar CrossRef Search ADS © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Error analysis of a variational multiscale stabilization for convection-dominated diffusion equations in two dimensions
JF - IMA Journal of Numerical Analysis
DO - 10.1093/imanum/drx027
DA - 2017-06-19
UR - https://www.deepdyve.com/lp/oxford-university-press/error-analysis-of-a-variational-multiscale-stabilization-for-j3e7EV7ZOx
SP - 1
EP - 1253
VL - Advance Article
IS - 3
DP - DeepDyve