
==== Front
BIT Numer Math
BIT Numer Math
Bit. Numerical Mathematics
0006-3835
1572-9125
Springer Netherlands Dordrecht

1035
10.1007/s10543-024-01035-8
Research Paper
Super-localized orthogonal decomposition for convection-dominated diffusion problems
Bonizzoni Francesca 1
http://orcid.org/0000-0002-9838-6321
Freese Philip philip.freese@tuhh.de

2
Peterseim Daniel 3
1 https://ror.org/01nffqt88 grid.4643.5 0000 0004 1937 0327 MOX-Dipartimento di Matematica, Politecnico di Milano, Piazza Leonardo da Vinci 32, 20133 Milan, Italy
2 grid.6884.2 0000 0004 0549 1777 Institute of Mathematics, Hamburg University of Technology, Am Schwarzenberg-Campus 3, 21073 Hamburg, Germany
3 https://ror.org/03p14d497 grid.7307.3 0000 0001 2108 9006 Institute of Mathematics & Centre for Advanced Analytics and Predictive Sciences (CAAPS), University of Augsburg, Universitätsstr. 12a, 86159 Augsburg, Germany
5 8 2024
5 8 2024
2024
64 3 333 12 2023
29 7 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
This paper presents a novel multi-scale method for convection-dominated diffusion problems in the regime of large Péclet numbers. The method involves applying the solution operator to piecewise constant right-hand sides on an arbitrary coarse mesh, which defines a finite-dimensional coarse ansatz space with favorable approximation properties. For some relevant error measures, including the L2-norm, the Galerkin projection onto this generalized finite element space even yields ε-independent error bounds, ε being the singular perturbation parameter. By constructing an approximate local basis, the approach becomes a novel multi-scale method in the spirit of the Super-Localized Orthogonal Decomposition (SLOD). The error caused by basis localization can be estimated in an a posteriori way. In contrast to existing multi-scale methods, numerical experiments indicate ε-robust convergence without pre-asymptotic effects even in the under-resolved regime of large mesh Péclet numbers.

Keywords

Convection-dominated diffusion
Numerical homogenization
Multi-scale method
Super-localization
Singularly perturbed
Mathematics Subject Classification

65N12
65N15
65N30
35B25
http://dx.doi.org/10.13039/100010663 H2020 European Research Council 3597 Bonizzoni Francesca Daniel Peterseim865751 Freese Philip Technische Universität Hamburg (3140)Open Access funding enabled and organized by Projekt DEAL.

issue-copyright-statement© BIT Foundation 2024
==== Body
pmcIntroduction

This paper studies the numerical solution of the following singularly perturbed convection-diffusion problem in a bounded, simply connected polygonal domain Ω⊂Rd with dimension d=1,2,3. Given some small diffusivity 0<ε≪1, a bounded, divergence-free and curl-free velocity field b as well as an external force f, we look for u such that the boundary value problem1.1 -εΔu+b·∇u=finΩu=0on∂Ω

holds in suitably weak sense.

This fairly simple model problem appears to be very challenging for classical Galerkin finite element methods (FEMs) and related schemes when the ratio of the convection rate over the diffusion is large, that is, for a large Péclet number Pe=bL∞(Ω)ε-1. In this regime, the solution u typically develops exponential and parabolic layers at the boundary (and possibly interior layers in the presence of inhomogeneous Dirichlet data). Unless the width h of the FE mesh resolves the characteristic length scale 1/Pe≈ε of these layers, FE approximations show spurious oscillations. To avoid this unstable pre-asymptotic behavior, a minimal resolution condition of the form hPe≲1 is typically required. However, in many relevant practical applications, ε may be so small that such conditions are unfeasible.

The circumvention or at least relaxation of this resolution condition has been the subject of intensive research in the past few decades. We refer to the monograph [62] for a detailed overview of the subject. Several branches of solution strategies have been developed. One is based on mesh refinement or grading toward layers [3, 31, 54, 55]. The more popular alternative, in particular in engineering communities, is the class of stabilized methods. Roughly speaking, these approaches change the model on the continuum or discrete level by adding artificial diffusion along the negative velocity field (upwinding). Among the extensive number of existing approaches in this context, we mention the streamline upwind/Petrov–Galerkin method [19] (also known as streamline diffusion method - see, e.g., [45]), the Galerkin least-squares method [41], the Douglas–Wang Galerkin method [34], algebraic flux correction [5, 48], discontinuous Petrov–Galerkin methods [26, 51], hybridizable discontinuous Galerkin methods [61], residual-free bubble methods [18, 21], nonconforming stabilized virtual element methods [7] and edge-based methods with additional nonlinear diffusion [4].

It has been observed that many of these stabilized schemes are strongly related to multi-scale methods, marking a third class of approaches to tackle strong convection [43]. The essential idea of multi-scale methods is to resolve the fine-scale features such as strong gradients in the layers by locally precomputed generalized FE shape functions. Prime examples are variational multi-scale methods (VMS) [42, 44, 49, 56], multi-scale FEMs [25, 58], multi-scale hybrid-mixed methods [38], multi-scale discontinuous Galerkin methods [23, 46], multi-scale virtual element methods [64], multi-scale stabilization methods [20, 22], stabilization procedures by means of sub-grid scale [24], energy minimizing generalized multi-scale methods [65], or the multi-scale method for time-dependent convection-dominated problems recently proposed in [63].

Although many of the approaches mentioned so far have been empirically successful in applications and certainly improved on the stability of standard FEMs, ε-independent behaviour is hardly observed for large mesh Péclet numbers hPe≫1. This statement also applies to the Localized Orthogonal Decomposition (LOD) method which originated from VMS and is often referred to as numerical homogenization (for an overview on the topic, see [2, 53, 57]). On an ideal level, the methodology realizes a prescribed projection of the unknown solution onto a discrete space (other than the Ritz projection) and hence allows best approximation results in suitable norms independent of the Péclet number. However, existing practical versions based on the localization of the fine-scale Green’s function [27, 43, 60] do suffer from strong convection. Although for moderate mesh Péclet numbers exponential decay results of [47, 52] for the fine-scale Green’s function still apply, they deteriorate with increasing mesh Péclet numbers as outlined in [50]. This prevents the construction of a localized basis by means of fine-scale correctors and limits the practical relevance of the approach.

An alternative localization strategy was recently proposed in [40] for the pure diffusion problem and then extended to indefinite and non-Hermitian problems in [36]. As described in [2], the LOD (and also the VMS) implicitly computes its problem-adapted ansatz space by applying the solution operator to some classical FE spaces on coarse meshes. For the specific choice of piecewise constants, the coarse space is simply given by the span of functions A-11T, A-1 denoting the solution operator and 1T being the characteristic function of the element T ranging into a coarse mesh TH. We refer to the Galerkin projection method on such ansatz space as ideal method. The novel localization strategy aims to identify local linear combinations of characteristic functions in such a way that the spread of the response under the solution operator is minimized. Since for the diffusion model problem this strategy yields a super-exponentially decaying localization error (as compared to the exponentially decaying localization error in classical LOD) the resulting practical method is referred to as Super-Localized Orthogonal Decomposition (SLOD).

The present paper shows that the super-localization strategy is not merely an amplification of the fine-scale Green’s function, but allows localization in applications where it has not been observed before. We generalize the SLOD methodology to convection-diffusion problems with large Péclet number. The SLOD approximation error comprises two contributions: the discretization error of the ideal method and the localization error. As such, the error analysis consists of two major steps. The key result to bound the first contribution is contained in Lemma 2.1, where a priori estimates for the continuous convection-diffusion problem with divergence-free and curl-free velocity field are proved. Thanks to this result, ε-explicit (and in particular cases, even ε-independent) error upper bounds for the ideal method are derived. The second contribution, instead, is proved to be proportional to the computable quantity σ (5.5), which reflects the worst-case localization error.

Notably, the SLOD basis functions display an ε-robust behaviour in numerical experiments. Indeed, as ε gets smaller, they seem not to be affected by oscillations nor an increase of their support seems necessary (see Figs. 2 and 3 for a representation in the one- and two-dimensional frameworks). This represents a major improvement with respect to both classical LOD and the state-of-the-art multi-scale method in [50]. From a practical point of view, this translates into significant computational savings, which in turn makes computations possible even in the three-dimensional framework (see Sect. 7.2 for 3D numerical experiments).

The remainder of the paper is organized as follows. In Sect. 2 a detailed description of the problem of interest in its variational formulation is shown, and a priori upper bounds for the continuous solution of the convection-diffusion problem with divergence-free and curl-free velocity field are proven. An ideal numerical homogenization method based on the L2-orthogonal projection onto piecewise constants is introduced in Sect. 3. The core of the paper are Sects. 4 and 5, where the novel localization approach is presented, turned into a practically feasible method and the error analysis is carried out. Section 6 explains the SLOD algorithm and in Sect. 7 its performances are displayed by means of several two- and three-dimensional numerical experiments.

Model problem

Let Ω⊂Rd be a polygonal, simply connected domain with d=1,2,3 and let 0<ε≤1 be a singular perturbation parameter. Moreover, let b∈L∞(Ω;Rd) satisfy divb=0 and curlb=0. Let V:=H01(Ω) and define the bilinear form a:V×V→R by2.1 a(u,v):=ε∫Ω∇u·∇vdx+∫Ω(b·∇u)vdx

for all u,v∈V. Given some linear functional F∈V′:=H-1(Ω) on V then the weak formulation of the boundary value problem (1.1) seeks u∈V such that, for all v∈V,2.2 a(u,v)=F(v).

From now on, we assume that the right-hand side is a bit more regular than minimal, i.e., it is of the form F(∙):=f,∙L2(Ω) for some f∈L2(Ω). This additional regularity of the right-hand side will give rise to orders of approximations. We focus on the convection-dominated regime, namely, ε≪1 and Péclet number Pe=bL∞(Ω)ε-1≫1.

Remark 2.1

The method proposed below naturally applies to the case of non-constant diffusion coefficients, which may incorporate multi-scale features, i.e., the constant diffusivity ε can be replaced by a variable one of the form εA where A∈L∞(Ω;Rd×d) is symmetric and positive definite almost everywhere in Ω. Moreover, the method can be generalized to the case of convection-diffusion–reaction equations in a straight-forward way.

Since divb=0, integration by parts implies, for all v∈V,2.3 a(v,v)=ε∫Ω∇v2dx+∫Ω(b·∇v)vdx=εvV2,

where ∙V=∇∙L2(Ω) denotes the H1-seminorm, which is a norm in V. Moreover, for all u,v∈V, the application of Cauchy–Schwarz’s and Poincaré’s inequalities readily implies2.4 a(u,v)≤CauVvV,

for Ca=Ca(Ω,bL∞(Ω))=ε+CPbL∞(Ω)>0, where CP denotes the Poincaré constant. The Lax-Milgram theorem, the coercivity (2.3) and the boundedness (2.4) show that Problem (2.2) admits a unique solution u∈V that satisfies the ε-dependent stability estimate2.5 uV≤CaεFH-1(Ω).

For F(∙)=(f,∙)L2(Ω) and special velocities, the estimate can be sharpened. More importantly, in the weaker L2(Ω)-norm, even ε-independent stability results are possible. We refer to [28, Lemma 2.1] which covers the special case b=10⊤. The subsequent lemma generalizes [28, Lemma 2.1] to velocity fields fulfilling the following technical assumption:

Assumption 1

The velocity field b is both curl- and divergence-free. Moreover, there exists bmin>0 such that |b|≥bmin almost everywhere.

Assumption 1 seems rather restrictive but just extends the assumption on b such that we can find a potential field ψ with b=∇ψ. We introduce Cψ=expψL∞(Ω). This is crucial for Lemma 2.1 below. Often rigorous numerical analysis requires constant fields or implicitly assume that a result like in Lemma 2.1 below holds. In the numerical experiments, we show that our scheme also applies in the case of more general velocity fields b.

Our results are phrased in the ε-scaled norm of V2.6 ∙V,ε2:=ε∙V2+∙L2(Ω)2,

which is equivalent to the ∙V-norm for ε≤1, since for all v∈V there holds2.7 εvV≤vV,ε≤1+CP2vV.

Lemma 2.1

Let b satisfy Assumption 1, then the unique solution of (2.2) with f∈L2(Ω) satisfies the estimate:uV,ε≤Cψ(1-ε)bmin2+Cψ2(1-ε)bmin2fL2(Ω).

In particular, for ε≤12, there holdsuV,ε≤CstabfL2(Ω),

with Cstab positive and ε-independent.

Proof

Following Assumption 1 and using that Ω is simply connected, we deduce that b=∇ψ for some harmonic potential ψ. Now, consider the transformed dependent variable v(x)=exp(-ψ(x))u(x), for all x∈Ω. To derive the strong formulation for v, we substitute u(x)=v(x)exp(ψ(x)) in (1.1). This yields:-εexp(ψ)Δv+2∇ψ·∇v+v∇ψ2+Δψ+exp(ψ)b·∇v+v∇ψ=f.

As the velocity field is curl- and divergence-free, we get ∇ψ=b and Δψ=0. Hence, after multiplying with vexp(-ψ) and using integration by parts as well as exploiting divb=0, we findεvV2+(1-ε)b2v,vL2(Ω)=f,exp(-ψ)vL2(Ω).

Consequently, we get the following estimateεvV2+(1-ε)bmin2vL2(Ω)2≤CψfL2(Ω)vL2(Ω)

with the constant Cψ=expψL∞(Ω). This yields a bound on the L2(Ω)-norm of v asvL2(Ω)≤Cψ(1-ε)bmin2fL2(Ω).

Eventually, to bound the L2(Ω)-norm of the solution u, we useuL2(Ω)=exp(ψ)vL2(Ω)≤CψvL2(Ω)≤Cψ2(1-ε)bmin2fL2(Ω)

The estimate on the ∙V-norm of the original solution u follows byεuV2=a(u,u)=f,uL2(Ω)≤fL2(Ω)uL2(Ω)≤CψfL2(Ω)vL2(Ω)≤Cψ2(1-ε)bmin2fL2(Ω)2

By combining the upper bounds on εuV2 and uL2(Ω), we derive the desired estimate. The ε-independent bound follows directly using 0<ε≤12 with Cstab=2bmin-2Cψbmin2+Cψ2.

□

Let us emphasize that the choice of the norm in Lemma 2.1 may be unusual or too weak in the context of convection dominated diffusion problems. However, we are not aware of any global estimates using stronger norms that include the directional derivative. Estimates using norms restricted to subdomains that avoid the layers may be possible, but are the subject of future research.

Remark 2.2

For the case of a convection-diffusion–reaction equation, the result from Lemma 2.1 is well known, but relies on the presence of the reaction term, see [62, Lemma 1.18]. In this case, as well as the special convection-diffusion case with b=10⊤ also (local) estimates on the directional derivative away from boundary layers are known, see [28, Lemma 1.2] and [62, Remark 1.19].

An ideal multi-scale method

This section introduces an ideal multi-scale method that identifies an approximation of the solution u in an operator-adapted ansatz space VH, whose construction is based on some (possibly coarse) FE mesh.

Let TH be a (triangular or quadrilateral) shape-regular quasi-uniform mesh (without hanging nodes) of the domain Ω, where H denotes the global mesh size of TH, namely, H=maxT∈THdiam(T). The degrees of freedom of the multi-scale method are associated with the mesh elements T∈TH via the characteristic functions 1T. Given the solution operator A-1:L2(Ω)→V that maps each right-hand-side function f∈L2(Ω) to the corresponding unique weak solution of problem (1.1) and the standard FE spaceP0(TH):=span1T|T∈TH

of TH-piecewise constants, the finite-dimensional subspace VH⊂V is given by3.1 VH:=A-1P0(TH)=spanA-11T|T∈TH.

Note that we could have chosen FE spaces other than P0(TH) for the approximation of the right-hand side. For example, the paper [27] considers discontinuous piecewise linears on simplicial meshes and [56] considers continuous piecewise linears with zero boundary condition. More generally, a finite-dimensional space of linear functionals on V could be considered. The authors in [50] implicitly use the Dirac delta functions δz for the interior vertices z of TH. Clearly this is only possible in one dimension and requires regularization in higher dimensions. While in two dimensions this was somewhat justifiable, the three-dimensional case seemed not to be tractable with this choice.

Let ΠH:L2(Ω)→P0(TH) denote the L2-orthogonal projection operator and note that, for all T∈TH, ΠHv|T is given byΠHv|T=1|T|∫Tvdx.

It is well-known that ΠH fulfills the following local stability and approximation properties (see [6, 59])3.2 ΠHvL2(T)≤vL2(T)for allv∈L2(T),

3.3 v-ΠHvL2(T)≤π-1H∇vL2(T)for allv∈H1(T).

Given the kernel W:=ker(ΠH|V) of ΠH when restricted to V, we may identify the decomposition3.4 V=VH⊕W,a(VH,W)=0,

which justifies the term orthogonal decomposition in the name of the method.

We will not use this decomposition explicitly in this paper, but our derivation is based on the composition of the L2 orthogonal projection ΠH and the solution operator A-1, which defines an ideal multi-scale method that maps right-hand sides f∈L2(Ω) into VH. Thus, uH=A-1ΠHf satisfies a(uH,v)=ΠHf,vL2(Ω) for all v∈V. Since uH is in VH by construction, it coincides with its Galerkin projection onto VH, i.e., uH∈VH is equivalently and uniquely characterized by the discrete variational problem3.5 a(uH,vH)=ΠHf,vHL2(Ω)

for all vH∈VH. Note that this is a non-standard projection of the solution u onto the discrete space VH. The method differs slightly from the classical Galerkin LOD or its Petrov–Galerkin variants. It equals the Galerkin projection and the abstract Petrov–Galerkin framework of [2] only for f∈P0(TH). For general f∈L2(Ω) it differs from the more established variants. In the pure diffusion case it equals the collocation variant discussed in [40].

In the following lemma we derive an ε-independent upper bound on the discretization error under Assumption 1, which is the motivation for the particular choice of the LOD variant (3.5).

Lemma 3.1

Let f∈Hs(Ω) with s∈[0,1], and b as in Assumption 1. Also assume that the assumption on ε from Lemma 2.1 is satisfied. Denote with u∈V and uH∈VH the unique solutions to (2.2) and (3.5), respectively. Then, there holds3.6 u-uHV,ε≤Cstabf-ΠHfL2(Ω)≤CCstabHsfHs(Ω),

where C, Cstab are ε- and H-independent positive constants, Cstab being introduced in Lemma 2.1.

Proof

Since u=A-1f and uH=A-1ΠHf we readily getu-uHV,ε=A-1f-A-1ΠHfV,ε=A-1(f-ΠHf)V,ε.

Lemma 2.1 provides an upper bound of the right-hand side. Altogether,u-uHV,ε≤Cstabf-ΠHfL2(Ω)≤CCstabHsfHs(Ω),

where the last inequality holds for all right-hand sides f∈Hs(Ω) with s∈[0,1]. □

Apart the exactness of the ideal method for f∈P0(TH), Lemma 3.1 above contains an error bound in the weaker L2(Ω)-norm that is independent of ε. First order convergence is predicted without a pre-asymptotic regime. The numerical experiments of the later sections will rather report second order and even ε-independent first order for the H1(Ω)-seminorm. A more abstract version of the estimate of (3.6) readsu-uHY≤‖A-1‖X→Yf-ΠHfX,

where ‖A-1‖X→Y refers to the norm of A-1 as a mapping between suitable Banach spaces X and Y. Choosing X=H-1(Ω) and Y=L2(Ω) or Y=H1(Ω) or Y=H1(ω) where ω⊂Ω excludes the boundary layers would pave the way to proving the numerically observed rates. However, we are not aware of any ε-independent bounds of the required operator norms.

Super-localization of basis functions

The canonical basis functions {A-11T|T∈TH} of the operator-adapted approximation space VH are non-local. To make the method practically feasible, localized basis functions have to be identified. The LOD provides a mechanism to construct an exponentially decaying basis that has been very successful in many applications. However, this is not the case when applied to convection-dominated problems, as we are interested here. More precisely, when applying the abstract theory of [2] the exponential decay property deteriorates as ε goes to 0, and the error estimate of error committed by computing a localized approximation of the exponentially decaying basis is only shown to behave like ε-1H-1-d/2exp(-cεℓ). This indicates that the localization parameter needs to grow algebraically in ε-1 to make this quantity small. This is in line with practical experience, documented, e.g., in [50]. Therein, the authors also discuss a possible improvement using anisotropic patches. However, the construction is based on point evaluation functionals and, hence, essentially limited to the one- and two-dimensional case.

This section presents an advanced localization strategy, which has superior properties, yielding, in particular, super-exponential decay of the localization error. The main idea stays in the identification of local TH-piecewise constant source terms that yield rapidly decaying (or even local) responses under the solution operator A-1 of the convection-dominated problem (1.1). This super-localization strategy, now known as the Super-Localized Orthogonal Decomposition (SLOD), has been first introduced in [40] for the second order elliptic partial differential equation -div(A∇u)=f, and subsequently extended to indefinite non-hermitian problems in [36].

For the subsequent derivation of the super-localization strategy, we need to introduce some notations. The local patch of level ℓ∈N of a union of elements S⊂Ω is given by:Nℓ(S):=⋃{T∈TH|T∩S≠∅}ℓ=1N1(Nℓ-1(S))ℓ=2,3,4,…

Let ℓ∈N be fixed, such that no patch coincides with the entire domain Ω. Given T∈TH, denoteω:=Nℓ(T) its ℓ-th order patch;

Vω:=v|ω|v∈V the restriction of V to the patch ω, equipped with the semi-norm ∙H1(ω) and the norm ∙H1(ω);

TH,ω:={K∈TH∩ω} the sub-mesh of TH with elements in ω;

ΠH,ω:L2(Ω)→P0(TH,ω) the L2-orthogonal projection onto P0(TH,ω).

Note that throughout the paper we will not distinguish between the functions in H01(ω) and their V-conforming extension by 0 to the entire domain Ω.

The (ideal) basis function φ=φT,ℓ,ε∈VH associated with the element T is given by the ansatzφ=A-1gwithg=gT,ℓ,ε:=∑K∈TH,ωcK1K,

for some coefficients (cK)K∈TH,ω that will be determined later. In particular, φ fulfils, for all v∈V,a(φ,v)=g,vL2(ω).

The Galerkin projection of φ onto the local subspace H01(ω) is the function φloc=φT,ℓ,εloc∈H01(ω) that satisfies, for all v∈H01(ω),4.1 aω(φloc,v)=g,vL2(ω),

where aω(·,·) denotes the restriction of the bilinear form a(·,·) to the subset ω. In general, the local function φloc is a poor approximation of the ideal function φ. Nevertheless, appropriate nontrivial choices of g (i.e., of coefficients (cK)K∈TH,ω) lead to highly accurate approximations in the energy norm.

The quantity we aim to minimize is the localization error φ-φloc, i.e. the error between the ideal basis function φ and its localized counterpart φloc. Note that the patch ω is a polytope. Hence, following [30, Theorem 31.31, Theorem 31.33], φloc∈H1+s(ω) for some s>12. The normal derivative is therefore integrable [29, Example 4.16], and we may derive the following lemma.

Lemma 4.1

(Variational characterization of the localization error) Let n be the outward normal of ω. There holds for all v∈Va(φ-φloc,v)=-ε∫∂ω\∂Ωn·∇φlocvds.

Moreover, we get a bound for the localization error as:φ-φlocV≤CCP2+1ℓHn·∇φlocL2(∂ω\∂Ω).

Proof

Let v∈V. Then, using the definitions of φ and φloc and integration by parts, there holds:a(φ-φloc,v)=a(φ,v)-a(φloc,v)=g,vL2(ω)-aω(φloc,v)=g,vL2(ω)-ε∫ω∇φloc·∇vdx-∫ω(b·∇φloc)vdx=g,vL2(ω)-∫ω-εΔφloc+(b·∇φloc)vdx-ε∫∂ω\∂Ω(n·∇φloc)vds=-ε∫∂ω\∂Ω(n·∇φloc)vds.

In particular, for eloc=φ-φloc, we findεelocV2=a(eloc,eloc)≤ε∫∂ω\∂Ω(n·∇φloc)elocds≤εn·∇φlocL2(∂ω\∂Ω)elocL2(∂ω\∂Ω).

Now, using the trace inequality (note that in the current setting we consider here, this dependence is of the form CℓH with some constant C>0 independent of H and ℓ), we getεelocV2≤εCℓHn·∇φlocL2(∂ω\∂Ω)elocH1(ω)≤εCℓHn·∇φlocL2(∂ω\∂Ω)elocH1(Ω).

Finally, from Poincarés inequality we findεelocV2≤εCCP2+1ℓHn·∇φlocL2(∂ω\∂Ω)elocV,

which yields the final result. □

Remark 4.1

In previous work [36, 40], the smallness of the normal derivative has been interpreted as the (almost) L2-orthogonality of g with respect to the space of convection-harmonic functions. Here, however, we directly use the smallness of the normal derivative, which makes the algorithm even simpler and avoids the sampling of the respective space of convection-harmonic functions.

Lemma 4.1 gives rise to the correct definition of the coefficients (cK)K∈TH,ω. Namely, we choose (cK)K∈TH,ω that minimize n·∇φlocL2(∂ω\∂Ω) subject to gL2(ω)=1. For the vector c consisting of the coefficients cK, this optimization task can be written in matrix form as follows:4.2 mincc⊤Ncsubject toc⊤Bc=1,

where the entries of the matrix N are (N)i,j=∫∂ω\∂Ω(n·∇φiloc)(n·∇φjloc)ds and B=HdI, I being the identity matrix. Equation (4.2) can be equivalently written as the following generalized eigenvalue problem: Find the smallest eigenvalue λ and the corresponding eigenvector c such that4.3 Nc=λBc.

More details on the precise choice of the basis functions is given in Sect. 6.

For the case of pure diffusion [40] there is strong numerical evidence that a quantity related to the L2(∂ω\∂Ω)-norm of the normal derivative decays super-exponentially in ℓ. Also in the presence of convection and even for high Péclet numbers, the numerical experiments in Fig. 1 show super-exponential decay of the smallest eigenvalue λ solution of problem (4.3), with respect to the localization parameter ℓ. As a consequence, there is a super-exponential decay of the normal derivative n·∇φlocL2(∂ω\∂Ω) of the corresponding localized basis functions in the localization parameter ℓ. For the limit ε→0, however, we believe that this decay deteriorates, and especially for the pure transport case, such a decay is not expected. However, in the presence of some diffusion in the model, the decay seems to be (super-)exponential. Henceforth, we assume that there exist g with gL2(ω)=1 and constants Csd(ε,H,ℓ)>0 depending on ε,H and ℓ, but independent of T, and C>0 independent of H,ℓ and T such that4.4 n·∇φlocL2(∂ω\∂Ω)≤Csd(ε,H,ℓ)exp-C(ε)ℓdd-1.

Fig. 1 All eigenvalues of the pair of matrices (N, B) for different values of the localization parameter ℓ, for a patch that does not reach the boundary . (left) Two-dimensional result for b=25+51+521⊤ and ε=2-11 on a coarse mesh H=2-4. (right) Three-dimensional result for b=27+51+5211⊤ and ε=2-6 on a coarse mesh with H=2-4

Remark 4.2

(SLOD basis in 1d) In the one-dimensional case, the boundary of the patches consists only of the two end points of the respective intervals, whereas we have three degrees of freedom for an order ℓ=1 patch. Thus, the minimization problem (4.2) can be solved exactly, which yields a vanishing normal derivative on both end points of the patches. Hence, from (4.4) for d=1, interpreting dd-1 as infinity, reveals a truly local basis function. This effect is also observed in Fig. 2, where we compare three different basis functions in VH for various values of ε and corresponding to the same mesh element T∈TH, namely A-11T (left); the basis function for L2-projection based LOD (center); the SLOD basis function φT,1,εloc (right).

Fig. 2 Solution to the convection-dominated problem with right-hand side 1T, i.e., A-11T (left); L2-projection based (global) LOD basis function (center); Novel SLOD basis function (right). Their corresponding L2-normalized right-hand sides are depicted in orange

Remark 4.3

(SLOD basis in 2d and 3d) While in the one-dimensional setting we were able to retrieve truly local basis functions, this is no longer true in higher dimensions. In Fig. 3 we depict the basis functions φT,4,εloc for an element T whose patch does not reach the global boundary for ε=2-9 and ε=2-11. The velocity field b is given as b=25+51+521⊤. Moreover, the figure shows the response of the solution operator to the indicator function 1T that corresponds to T. It is clearly visible, that the SLOD basis functions decay very fast, especially in comparison to the ideal basis functions of the space (3.1).

Fig. 3 Absolute value of SLOD basis (left) and solution of A-11T (right) on 4-th order (interior) patch, for b=25+51+521⊤ and different values of ε

Super-localized multi-scale method and error analysis

Within this section we turn the method (3.5) based on the ideal operator-adapted ansatz subspace VH⊂V into a feasible numerical scheme, by means of the super-localization strategy introduced above.

Let the localization parameter ℓ be fixed. We define the ansatz space of the super-localized method as the span of the SLOD basis functions φT,ℓ,εloc as T varies in the coarse grid TH, namely:5.1 VHℓ:=φT,ℓ,εloc|T∈TH⊂V.

The approximate solution provided by the SLOD method is the Galerkin projection in the space VHℓ of the convection-dominated problem at hand with perturbed right-hand side ΠHf. In particular, the SLOD approximation to (2.2) is the function uHℓ∈VHℓ such that, for all vHℓ∈VHℓ,5.2 a(uHℓ,vHℓ)=ΠHf,vHℓL2(Ω).

We emphasize that we can expand ΠHf∈P0(TH) in the basis of right-hand sides gT,ℓ,ε|T∈TH, which yields5.3 ΠHf=∑T∈THcTgT,ℓ,ε.

A minimal requirement for the stability and convergence of the Galerkin method (5.2) is that the set of functions gT,ℓ,ε|T∈TH spans P0(TH) in a stable way. Numerically, this is ensured as described in Sect. 6. For the subsequent theoretical analysis, we make the following assumption.

Assumption 2

The set gT,ℓ,ε|T∈TH is a Riesz basis of P0(TH), i.e., there exists a constant Crb(ε,H,ℓ), depending only polynomially on H and ℓ, such that, for all {cT}T∈TH⊂R, there holds5.4 Crb-1(ε,H,ℓ)∑T∈TH|cT|2≤∑T∈THcTgT,ℓ,εL2(Ω)2≤Crb(ε,H,ℓ)∑T∈TH|cT|2.

In the following theorem we derive an a priori error estimate for the solution to problem (5.2). The upper bound is explicit in the quantity5.5 σ(ε,H,ℓ):=maxT∈THn·∇φT,ℓ,εlocL2(∂ω\∂Ω)

which reflects the worst-case localization error.

Remark 5.1

(Exponential decay of classical LOD) For moderate mesh Péclet number, the quantity σ(ε,H,ℓ) in (5.5) decays exponentially in the localization parameter ℓ (see [40, Appendix A] for the proof in the pure diffusion case). In particular, one can recover the a priori error estimate with rates similar to those for the LOD theory as in [2, 27, 50].

Theorem 5.1

(Convergence of the SLOD method) Let Assumption 2 and Assumption 1 be satisfied. Then, there exists a constant C>0 independent of H,ℓ,ε such that, for all f∈Hs(Ω) with s∈[0,1], there holds5.6 u-uHℓV,ε≤CCstabf-ΠHfL2(Ω)+σ(ε,H,ℓ)εℓHCrb(ε,H,ℓ)1/2ℓd/2fL2(Ω)≤CHsfHs(Ω)+σ(ε,H,ℓ)εℓHCrb(ε,H,ℓ)1/2ℓd/2fL2(Ω),

where Crb(ε,H,ℓ) and σ(ε,H,ℓ) are defined in Assumption 2 and (5.5), respectively.

Proof

By triangular inequality, we get:5.7 u-uHℓV,ε≤u-uHV,ε+uH-uHℓV,ε.

The first term in (5.7) represents the discretization error of the ideal multi-scale method, and its upper bound is given by Lemma 3.1. We consider now the second term in (5.7), which represents the localization error. Observe that uH solves the continuous equation for right-hand side ΠHf. As a consequence, the SLOD solution uHℓ is the Galerkin approximation of uH in the finite dimensional space VHℓ. Using the norm equivalence (2.7) and applying Céa’s Lemma, we getuH-uHℓV,ε≲uH-uHℓV≲1εinfvHℓ∈VHℓuH-vHℓV,

where the notation x≲y means x≤cy with c positive constant independent of the mesh size parameter H, the localization parameter ℓ and the diffusion coefficient ε. Given the expansion of ΠHf in the basis gT,ℓ,ε|T∈TH, namely, ΠHf=∑T∈THcTgT,ℓ,ε, we can express uH asuH=∑T∈THcTA-1gT,ℓ,ε=∑T∈THcTφT,ℓ,ε.

For the particular choice vHℓ=∑T∈THcTφT,ℓ,εloc, we obtain that e:=uH-vHℓ∈V fulfils:eV2=1εa(uH-vHℓ,e)=1ε∑T∈THcTa(φT,ℓ,ε-φT,ℓ,εloc,e)=-∑T∈THcT∫∂ωn·∇φT,ℓ,εloceds≤∑T∈THcTn·∇φT,ℓ,εlocL2(∂ω)CℓHeH1(ω)≤σ(ε,H,ℓ)CℓH∑T∈THcTeH1(ω),

where we employed Lemma 4.1 in the third equality and (5.5) in the last inequality. For simplicity, we omit the dependence of σ and Crb on ε,H and ℓ in the rest of the proof. As a consequence, thanks to Assumption 2, (5.3), the Poincaré inequality and (3.2), there holds:eV2≲σℓH∑T∈THcTeH1(ω)≲σℓH∑T∈THcT2∑T∈THeH1(ω)2≲σℓHCrb1/2ΠHfL2(Ω)Colℓd/2eH1(Ω)≲σℓHCrb1/2fL2(Ω)Colℓd/2eV,

where Col2ℓd bounds the number of patches containing a fixed mesh element. In particular, we have proved thateV≲σℓHCrb1/2ℓd/2fL2(Ω),

so that the estimate (5.6) follows. □

As previously observed for the ideal multi-scale method, upper bounds on the SLOD error could be derived in the abstract setting A-1:X→Y, for suitable Banach spaces X and Y.

We point out that in the case of a piecewise constant right-hand side f, the first term in (5.6) vanishes. From our assumption on the decay of the normal derivative in (4.4), we deduce that the ε-dependence of the second expression is dominated by the exponentially decaying quantity σ(ε,H,ℓ). Moreover, we derive that the localization condition ℓ≳log(εH)d-1d guarantees convergence of the SLOD with order H. In the limit ε→0, this would still result in a dense matrix.

Numerical implementation and stable selection of basis

This section discusses the implementation of the proposed numerical method, with particular attention to the computation of a basis {φT,ℓ,εloc|T∈TH} for the ansatz space VHℓ that is associated with a basis {gT,ℓ,εloc|T∈TH} of P0(TH) via  (4.1). The Riesz stability of the basis in the sense of Assumption 2 must be respected. However, there is still no a priori guarantee that this will hold, but we believe that the approach below helps to get a Riesz basis, which may be checked a posteriori.

For simplicity, we take Ω as the unit hypercube in d dimensions, i.e., Ω=(0,1)d, discretized by means of a quadrilateral mesh TH. Given ℓ≥1, we choose an element T∈TH and consider the corresponding patch ω=Nℓ(T). In a first step, for each element K∈TH,ω in the patch mesh, we compute the response of the solution operator restricted to the patch, denoted by Aω-1, to its characteristic function 1K, i.e., Aω-11K. By construction, the target basis function φT,ℓ,εloc is in the span of these #TH,ω≈ℓd local responses. In a second step, we search for the function φT,ℓ,εloc∈spanAω-11K|K∈TH,ω in this low-dimensional space by minimizing normal derivatives subject to a unit mass constraint as presented in (4.2). The corresponding eigenvector (cK)K∈TH,ω contains the coefficients of the expansion of φT,ℓ,εloc in terms of the local responses. At the same time, the coefficients are the values of gT,ℓ,εloc in the elements of the patch.

Unfortunately, the smallest eigenvalue may not be simple, or there might be a cluster of small eigenvalues. Then a particular choice of eigenfunction may not always be favourable with regard to the global stability of the basis in the sense of Assumption 2. Especially for patches that touch the boundary of the global domain Ω [40, Appendix B] and in the regime of high convection, an additional optimization step ensures linear independence of the functions computed in different patches. For this purpose, we incorporate eigenfunctions associated with a certain range of the lowermost eigenvalues. Given all eigenvalues λ1≤λ2≤⋯≤λ#TH,ω and some parameter p≥1, we choose all indices 1≤i≤#TH,ω so that6.1 λiλ#TH,ω≤maxλ1λ#TH,ω1p,1e-10,

and we denote the resulting set of indices by I. The choice p=1 reflects the case where only the smallest (potentially multiple) eigenvalue is used, and thus we use p>1 in our implementation.

Among these candidate functions with close to minimal normal derivative at the boundary of the patch, we choose the one that maximizes a weighted L2(ω)-norm under the unit mass constraint. The piecewise constant weight function is zero in the central element T and grows in a b-dependent way with a certain distance from the central element. Let us introduce the midpoints mT,mK∈Rd of the central element and an element of the patch, respectively. We define the relation between elements asrel(T,K):=H-1(mK-mT)∈Zd,

and introduce for each element K∈TH,ω\T its weight by6.2 wK:=rel(T,K)-b(mT)b(mT)2ℓ∞w.

Here pw≥1 is a parameter that needs to be chosen. With (6.2) we ensure that the elements in the b direction are penalized less, which accounts for the natural shift of the basis due to convection. For a realization of an order 1 patch, see Fig. 4.Fig. 4 Weights wK for an element that does not reach the boundary, for constant velocity b as given in (7.1), and pw=2

Eventually, we search the function in the space of the previously selected candidate right-hand sides spangT,ℓ,ε,i|i∈I that minimizes a weighted L2(ω)-norm subject to the unit mass constraint. This constraint minimization is realized by computing the smallest eigenvalue of the symmetric positive definite matrix6.3 1gT,ℓ,ε,iL2(ω)gT,ℓ,ε,jL2(ω)∑K∈TH,ω∫KwKgT,ℓ,ε,igT,ℓ,ε,jdxi,j∈I.

In this way, we compute for every element T of the coarse mesh TH the basis function φT,ℓ,εloc and hence build the space VH,ℓ. From our numerical experiments, the choices p=3 and pw=2 produce good results. In Algorithm 1 we detail the full algorithm for the computation of the super-localized basis.

Algorithm 1 Basis selection

Although our approach gives good results in the numerical experiments in Sect. 7, we believe that it can be improved. In particular, the choice of the two-step optimization leaves room for improvement to ensure favorable locality and stability at the same time. Nevertheless, this is the best optimization we have found so far. Another way to achieve effective stabilization may be to use partiton of unity techniques, as suggested in [35].

Numerical experiments

In this section, we demonstrate the performance of our method. For this purpose, we briefly introduce the general configuration. All our experiments have been performed in Matlab. The computational domain Ω is given as a unit hypercube in d dimensions, i.e., Ω=(0,1)d. We introduce a fine quadrilateral mesh Th, which resolves the small parameter ε and is used to compute reference solutions and localized basis functions using the standard Galerkin FE method on the space of piecewise bilinear polynomials. In addition, we consider a coarse quadrilateral mesh TH as the target scale, which does not resolve ε. We emphasize that on this mesh, the degrees of freedom of the FEM correspond to the vertices of the mesh, while the degrees of freedom of the SLOD correspond to the elements. Thus, in principle, the SLOD has fewer degrees of freedom, but an increasing localization parameter results in a slightly denser matrix.

Two-dimensional experiment

We start by presenting two-dimensional experiments and compare our approach with the streamline upwind/Petrov–Galerkin (SUPG) method (see [34]) that we briefly recall below. Let UH denote the standard Galerkin FE space of piecewise bilinear polynomials on the coarse mesh TH. The SUPG approximation uHSUPG∈UH satisfies, for all vH∈UHBSUPG(uHSUPG,vH)=FSUPG(vH),

withBSUPG(uHSUPG,vH):=a(uHSUPG,vH)+δSUPG∑T∈THb·∇uHSUPG,b·∇vHL2(T)

andFSUPG(vH)=⟨f,vH⟩H-1(Ω)×H01(Ω)+δSUPG∑T∈THf,b·∇vHL2(T).

The symbol δSUPG denotes the stabilization parameter, and is chosen as δSUPG=H2b2.

Constant velocity

First, we follow the experiment from [50, Section 6]. We choose the right-hand side f≡1 and the constant velocity field b as7.1 b(x)=cos(0.7)sin(0.7)⊤.

The singular perturbation parameter ε is chosen to be 2-8. In this configuration, we expect boundary layers at the right and top boundaries. Moreover, in this situation the right-hand side is piecewise constant and hence, the first expression in our error estimate in (5.6) vanishes. Therefore, we observe the localization error.

The mesh size h of the fine mesh Th is chosen as h=2-10. Figure 5 shows the corresponding reference solution, its FE and SLOD approximation on a coarse mesh with H=2-4. The localization parameter ℓ is chosen equal to 1.Fig. 5 Reference solution computed on fine mesh with h=2-10 (a), FE approximation (b) and SLOD approximation with ℓ=1 (c) on coarse mesh with H=2-4 for the constant velocity field b given in (7.1), right-hand side f≡1 and ε=2-8

We observe that the SLOD resolves the layer, whereas the classic FE approximation suffers from severe instabilities. Figure 6 shows the convergence rates of the SLOD method for different localization parameters ℓ and coarse mesh sizes H as well as the error of the FE and SUPG methods.Fig. 6 Error in L2- (left) and ∙V-norm (right) for the constant velocity field b as in (7.1) and right-hand side f≡1 with ε=2-8

Fig. 7 Decay of the localization error in ∙V-norm versus the localization parameter ℓ for ε=2-8 and various values of H

The super-exponential convergence of (4.4) is numerically verified in Fig. 7. Unfortunately, our method shows inaccuracies for refinements in H. Most likely, these are due to the selection of the basis functions as discussed in Sect. 6 and hence an improvement in this selection process could lead to a more accurate method. However, we chose to ensure stability and possibly lose some accuracy in return.

Variable velocity and non-constant right-hand side

In this example, we consider a varying velocity field b, given as7.2 b(x)=21+45sin(4π(x1+1))⊤,for allx∈Ω.

It should be emphasized that the velocity field is divergence-free but not curl-free. Nevertheless, the numerical experiments show that the proposed multi-scale method can be satisfactorily applied beyond the assumptions of the theory. Since b varies in space, the penalization introduced in Sect. 6 varies with the different macroscopic cells in the coarse mesh. The flow in this example yields a boundary layer at the right boundary. We choose the non-constant right-hand side as f(x)=sin(πx1)cos(πx2). Thus, with respect to the error estimate in Theorem 5.1 the first expression in (5.6) does not vanish, and we expect a convergence of order one in the ∙V-norm as the right-hand side is regular enough. Figure 8 shows the reference solution as well as its FE and SLOD approximations. We again observe that the SLOD resolves the boundary layer, whereas the FEM delivers poor results.Fig. 8 Reference solution computed on fine mesh with h=2-10 (a), FE approximation (b) and SLOD approximation with ℓ=1 (c) on coarse mesh with H=2-4 for the variable velocity field b given in (7.2), right-hand side f=sin(πx1)cos(πx2) and ε=2-8

For non-constant right-hand sides, we expect improved approximation properties of the SLOD method (5.1) with ΠHf replaced by f. We will refer to such a method as SLOD-Galerkin. Note that the two methods produce different approximations and require different computational efforts. More specifically, since both methods search for an approximation in the same ansatz space VHℓ, they share the offline phase, namely the computation of the set of operator-adapted local basis functions. On the other hand, they differ in the online phase, where the actual approximation to the solution is computed. In particular, the SLOD method is more efficient online than the SLOD-Galerkin method.

The convergence in the L2(Ω) and ∙V norms for both the SLOD and the SLOD-Galerkin is shown in Fig. 9, where we observe second order convergence in the ∙V norm in H for both methods. In the L2(Ω) norm, the SLOD-Galerkin has a third order convergence, one order better than our proposed method, due to the extra effort to integrate the right-hand side more accurately.Fig. 9 Error in L2(Ω)- (left) and ∙V-norm (right) for the variable velocity field b as in (7.2) and right-hand side f=sin(πx1)cos(πx2). We consider the parameter ε=2-8 and the proposed SLOD method as well as the SLOD-Galerkin (SLOD-G)

Three-dimensional experiment

As mentioned before, the variational multi-scale stabilization method from [50] only works in one or two dimensions. Here we show that the super-localized variant is capable to approximate the solution even in a three-dimensional setup. In this configuration we again choose a constant velocity field b, which is given as7.3 b=π4111⊤.

The constant right-hand side is f≡1. We expect a boundary layer around the top right corner at 111⊤. In the three-dimensional setting we compute the reference solution on a mesh with h=2-6, which resolves the chosen ε=2-5. Figure 10 shows the L2(Ω)- and ∙V-norm errors obtained for the SLOD method. As in our first experiment, due to the constant right-hand side we observe the localization error.Fig. 10 Error in the L2(Ω)- (left) and ∙V-norm (right) for the constant velocity field b as in (7.3), right-hand side f≡1 and ε=2-5, in the three-dimensional case

From (4.4), for d=3, we deduce that the localization error behaves like exp(-Cℓ1.5). Figure 11 illustrates this super-exponential decay in the ∙V-norm.Fig. 11 Super-exponential decay of the ∙V-norm in localization parameter ℓ for the three-dimensional experiment

Concluding remarks and future developments

We have presented a novel multi-scale method for convection-dominated problems. The method follows the LOD framework and employs a novel super-localization strategy. The resulting SLOD significantly improves previous attempts to tackle convection-dominated problems in the under-resolved regime of large mesh Péclet numbers. While previously the SLOD largely improved the performance of already very efficient methods for model diffusion and Helmholtz problems [36, 40], the present paper demonstrates the true potential of the super-localization idea to enlarge the class of problems tractable by multi-scale methods. The numerically observed ε-independent convergence in two as well as three-dimensional experiments is to some extent justified by numerical analysis involving a priori and a posteriori techniques.

Among the many promising future research directions are parameterized elliptic multi-scale problems to be treated by combining SLOD with model order reduction techniques, following the ideas in [1]. It has been shown in [9] that this is conceptually possible for problems involving convection. An interesting and relevant application from the physical point of view are wave propagation and scattering problems in highly heterogeneous structures. For such an objective, the SLOD has been recently proposed in [36], and model order reduction techniques for the parametric-in-frequency problem have recently been presented (see, e.g., [13–17]). Moreover, the possible improvement of numerical stochastic homogenization methods [32, 33, 37, 39] and uncertainty quantification techniques [8, 10–12] will be analyzed.

Acknowledgements

The authors gratefully acknowledge Gabriel Barrenechea for fruitful discussion on the stability properties of the convection-dominated boundary value problem. We also thank the anonymous reviewers for their constructive comments, which helped us to improve the manuscript substantially. In particular, on the detailed suggestion of one reviewer, the result in the presentation of Lemma 2.1 has been clarified and significantly extended.

Funding

Open Access funding enabled and organized by Projekt DEAL.

Declarations

Conflict of interest

The authors declare no Conflict of interest.

The work of all authors is part of a project that has received funding from the European Research Council ERC under the European Union’s Horizon 2020 research and innovation program (Grant agreement No. 865751).

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Abdulle A Henning P A reduced basis localized orthogonal decomposition J. Comput. Phys. 2015 295 379 401
Abdulle, A., Henning, P.: A reduced basis localized orthogonal decomposition. J. Comput. Phys. 295, 379–401 (2015)
2. Altmann R Henning P Peterseim D Numerical homogenization beyond scale separation Acta Numer. 2021 30 1 86
Altmann, R., Henning, P., Peterseim, D.: Numerical homogenization beyond scale separation. Acta Numer. 30, 1–86 (2021)
3. Bakhvalov NS On the optimization of the methods for solving boundary value problems in the presence of a boundary layer Zh. Vychisl. Mat. i Mat. Fiz. 1969 9 4 841 859
Bakhvalov, N.S.: On the optimization of the methods for solving boundary value problems in the presence of a boundary layer. Zh. Vychisl. Mat. i Mat. Fiz. 9(4), 841–859 (1969)
4. Barrenechea GR Burman E Karakatsani F Edge-based nonlinear diffusion for finite element approximations of convection-diffusion equations and its relation to algebraic flux-correction schemes Numer. Math. 2017 135 2 521 545 28615743
Barrenechea, G.R., Burman, E., Karakatsani, F.: Edge-based nonlinear diffusion for finite element approximations of convection-diffusion equations and its relation to algebraic flux-correction schemes. Numer. Math. 135(2), 521–545 (2017)28615743
5. Barrenechea GR John V Knobloch P Analysis of algebraic flux correction schemes SIAM J. Numer. Anal. 2016 54 4 2427 2451
Barrenechea, G.R., John, V., Knobloch, P.: Analysis of algebraic flux correction schemes. SIAM J. Numer. Anal. 54(4), 2427–2451 (2016)
6. Bebendorf M A note on the Poincaré inequality for convex domains Zeitschrift für Anal. und ihre Anwendungen 2003 22 4 751 756
Bebendorf, M.: A note on the Poincaré inequality for convex domains. Zeitschrift für Anal. und ihre Anwendungen 22(4), 751–756 (2003)
7. Berrone S Borio A Manzini G SUPG stabilization for the nonconforming virtual element method for advection-diffusion-reaction equations Comput. Methods Appl. Mech. Eng. 2018 340 500 529
Berrone, S., Borio, A., Manzini, G.: SUPG stabilization for the nonconforming virtual element method for advection-diffusion-reaction equations. Comput. Methods Appl. Mech. Eng. 340, 500–529 (2018)
8. Bonizzoni F Buffa A Nobile F Moment equations for the mixed formulation of the Hodge Laplacian with stochastic loading term IMA J. Numer. Anal. 2013 34 4 1328 1360
Bonizzoni, F., Buffa, A., Nobile, F.: Moment equations for the mixed formulation of the Hodge Laplacian with stochastic loading term. IMA J. Numer. Anal. 34(4), 1328–1360 (2013)
9. Bonizzoni F Hauck M Peterseim D A reduced basis super-localized orthogonal decomposition for reaction-convection-diffusion problems J. Comput. Phys. 2024 499 112698
Bonizzoni, F., Hauck, M., Peterseim, D.: A reduced basis super-localized orthogonal decomposition for reaction-convection-diffusion problems. J. Comput. Phys. 499, 112698 (2024)
10. Bonizzoni, F., Nobile, F.: Perturbation analysis for the stochastic Darcy problem. In: ECCOMAS 2012-European Congress on Computational Methods in Applied Sciences and Engineering, pp. 3926–3933 (2012)
11. Bonizzoni F Nobile F Perturbation analysis for the Darcy problem with log-normal permeability SIAM/ASA J. Uncertain. Quantif. 2014 2 1 223 244
Bonizzoni, F., Nobile, F.: Perturbation analysis for the Darcy problem with log-normal permeability. SIAM/ASA J. Uncertain. Quantif. 2(1), 223–244 (2014)
12. Bonizzoni F Nobile F Regularity and sparse approximation of the recursive first moment equations for the lognormal Darcy problem Comput. Math. Appl. 2020 80 12 2925 2947
Bonizzoni, F., Nobile, F.: Regularity and sparse approximation of the recursive first moment equations for the lognormal Darcy problem. Comput. Math. Appl. 80(12), 2925–2947 (2020)
13. Bonizzoni F Nobile F Perugia I Convergence analysis of Padé approximations for Helmholtz frequency response problems ESAIM Math. Modell. Numer. Anal. 2018 52 4 1261 1284
Bonizzoni, F., Nobile, F., Perugia, I.: Convergence analysis of Padé approximations for Helmholtz frequency response problems. ESAIM Math. Modell. Numer. Anal. 52(4), 1261–1284 (2018)
14. Bonizzoni F Nobile F Perugia I Pradovera D Fast least-Squares Padé approximation of problems with normal operators and meromorphic structure Math. Comput. 2020 89 323 1229 1257
Bonizzoni, F., Nobile, F., Perugia, I., Pradovera, D.: Fast least-Squares Padé approximation of problems with normal operators and meromorphic structure. Math. Comput. 89(323), 1229–1257 (2020)
15. Bonizzoni F Nobile F Perugia I Pradovera D Least-Squares Padé approximation of parametric and stochastic Helmholtz maps Adv. Comput. Math. 2020 46 3 1 28
Bonizzoni, F., Nobile, F., Perugia, I., Pradovera, D.: Least-Squares Padé approximation of parametric and stochastic Helmholtz maps. Adv. Comput. Math. 46(3), 1–28 (2020)
16. Bonizzoni, F., Pradovera, D.: Shape optimization for a noise reduction problem by non-intrusive parametric reduced modeling. In: 14th WCCM-ECCOMAS Congress 2020, vol. 700 (2021)
17. Bonizzoni F Pradovera D Ruggeri M Rational-based model order reduction of Helmholtz frequency response problems with adaptive finite elementsnapshots Math. Eng. 2023 5 4 1 38
Bonizzoni, F., Pradovera, D., Ruggeri, M.: Rational-based model order reduction of Helmholtz frequency response problems with adaptive finite elementsnapshots. Math. Eng. 5(4), 1–38 (2023)
18. Brezzi F Marini D Süli E Residual-free bubbles for advection-diffusion problems: the general error analysis Numer. Math. 2000 85 1 31 47
Brezzi, F., Marini, D., Süli, E.: Residual-free bubbles for advection-diffusion problems: the general error analysis. Numer. Math. 85(1), 31–47 (2000)
19. Brooks AN Hughes TJ Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations Comput. Methods Appl. Mech. Eng. 1982 32 1 199 259
Brooks, A.N., Hughes, T.J.: Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Eng. 32(1), 199–259 (1982)
20. Calo VM Chung ET Efendiev Y Leung WT Multiscale stabilization for convection-dominated diffusion in heterogeneous media Comput. Methods Appl. Mech. Eng. 2016 304 359 377
Calo, V.M., Chung, E.T., Efendiev, Y., Leung, W.T.: Multiscale stabilization for convection-dominated diffusion in heterogeneous media. Comput. Methods Appl. Mech. Eng. 304, 359–377 (2016)
21. Cangiani A Süli E Enhanced residual-free bubble method for convection-diffusion problems Int. J. Numer. Meth. Fluids 2005 47 10–11 1307 1313
Cangiani, A., Süli, E.: Enhanced residual-free bubble method for convection-diffusion problems. Int. J. Numer. Meth. Fluids 47(10–11), 1307–1313 (2005)
22. Chung ET Efendiev Y Leung WT Multiscale stabilization for convection-diffusion equations with heterogeneous velocity and diffusion coefficients Comput. Math. Appl. 2020 79 8 2336 2349
Chung, E.T., Efendiev, Y., Leung, W.T.: Multiscale stabilization for convection-diffusion equations with heterogeneous velocity and diffusion coefficients. Comput. Math. Appl. 79(8), 2336–2349 (2020)
23. Chung ET Leung WT A sub-grid structure enhanced discontinuous Galerkin method for multiscale diffusion and convection-diffusion problems Commun. Comput. Phys. 2013 14 2 370 392
Chung, E.T., Leung, W.T.: A sub-grid structure enhanced discontinuous Galerkin method for multiscale diffusion and convection-diffusion problems. Commun. Comput. Phys. 14(2), 370–392 (2013)
24. Codina R Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods Comput. Methods Appl. Mech. Eng. 2000 190 13 1579 1599
Codina, R.: Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods. Comput. Methods Appl. Mech. Eng. 190(13), 1579–1599 (2000)
25. Degond P Lozinski A Muljadi BP Narski J Crouzeix-Raviart MsFEM with bubble functions for diffusion and advection-diffusion in perforated media Commun. Comput. Phys. 2015 17 4 887 907
Degond, P., Lozinski, A., Muljadi, B.P., Narski, J.: Crouzeix-Raviart MsFEM with bubble functions for diffusion and advection-diffusion in perforated media. Commun. Comput. Phys. 17(4), 887–907 (2015)
26. Demkowicz L Gopalakrishnan J Niemi AH A class of discontinuous Petrov-Galerkin methods. Part iii: adaptivity Appl. Numer. Math. 2012 62 4 396 427
Demkowicz, L., Gopalakrishnan, J., Niemi, A.H.: A class of discontinuous Petrov-Galerkin methods. Part iii: adaptivity. Appl. Numer. Math. 62(4), 396–427 (2012)
27. Elfverson, D.: A discontinuous Galerkin multiscale method for convection-diffusion problems. arXiv preprint arXiv:1509.03523 (2015)
28. Eriksson K Johnson C Adaptive streamline diffusion finite element methods for stationary convection-diffusion problems Math. Comp. 1993 60 201 167 188
Eriksson, K., Johnson, C.: Adaptive streamline diffusion finite element methods for stationary convection-diffusion problems. Math. Comp. 60(201), 167–188 (1993)
29. Ern A Guermond JL Finite elements I-Approximation and interpolation Texts in Applied Mathematics 2021 Cham Springer 45 89
Ern, A., Guermond, J.L.: Finite elements I-Approximation and interpolation. In: Texts in Applied Mathematics, vol. 72, pp. 45–89. Springer, Cham (2021). 10.1007/978-3-030-56341-7
30. Ern A Guermond JL Finite elements II–Galerkin approximation, elliptic and mixed PDEs, Texts in Applied Mathematics 2021 Cham Springer
Ern, A., Guermond, J.L.: Finite elements II–Galerkin approximation, elliptic and mixed PDEs, Texts in Applied Mathematics, vol. 73. Springer, Cham (2021). 10.1007/978-3-030-56923-5
31. Farrell PA Hegarty AF Miller JJH O’Riordan E Shishkin GI Robust computational techniques for boundary layers, Applied Mathematics (Boca Raton) 2000 Boca Raton, FL Chapman & Hall/CRC
Farrell, P.A., Hegarty, A.F., Miller, J.J.H., O’Riordan, E., Shishkin, G.I.: Robust computational techniques for boundary layers, Applied Mathematics (Boca Raton), vol. 16. Chapman & Hall/CRC, Boca Raton, FL (2000)
32. Feischl M Peterseim D Sparse compression of expected solution operators SIAM J. Numer. Anal. 2020 58 6 3144 3164
Feischl, M., Peterseim, D.: Sparse compression of expected solution operators. SIAM J. Numer. Anal. 58(6), 3144–3164 (2020)
33. Fischer J Gallistl D Peterseim D A priori error analysis of a numerical stochastic homogenization method SIAM J. Numer. Anal. 2021 59 2 660 674
Fischer, J., Gallistl, D., Peterseim, D.: A priori error analysis of a numerical stochastic homogenization method. SIAM J. Numer. Anal. 59(2), 660–674 (2021)
34. Franca LP Frey SL Hughes TJ Stabilized finite element methods: I. Application to the advective-diffusive model Comput. Methods Appl. Mech. Eng. 1992 95 2 253 276
Franca, L.P., Frey, S.L., Hughes, T.J.: Stabilized finite element methods: I. Application to the advective-diffusive model. Comput. Methods Appl. Mech. Eng. 95(2), 253–276 (1992)
35. Freese P Hauck M Keil T Peterseim D A super-localized generalized finite element method Numer. Math. 2024 156 1 205 235
Freese, P., Hauck, M., Keil, T., Peterseim, D.: A super-localized generalized finite element method. Numer. Math. 156(1), 205–235 (2024)
36. Freese P Hauck M Peterseim D Super-localized orthogonal decomposition for high-frequency Helmholtz problems SIAM J. Sci. Comput. 2024 46 4 A2377 A2397
Freese, P., Hauck, M., Peterseim, D.: Super-localized orthogonal decomposition for high-frequency Helmholtz problems. SIAM J. Sci. Comput. 46(4), A2377–A2397 (2024)
37. Gallistl D Peterseim D Numerical stochastic homogenization by quasi-local effective diffusion tensors Commun. Math. Sci. 2019 17 3 637 651
Gallistl, D., Peterseim, D.: Numerical stochastic homogenization by quasi-local effective diffusion tensors. Commun. Math. Sci. 17(3), 637–651 (2019)
38. Harder C Paredes D Valentin F On a multiscale hybrid-mixed method for advective-reactive dominated problems with heterogeneous coefficients Multiscale Model. Simul. 2015 13 2 491 518
Harder, C., Paredes, D., Valentin, F.: On a multiscale hybrid-mixed method for advective-reactive dominated problems with heterogeneous coefficients. Multiscale Model. Simul. 13(2), 491–518 (2015)
39. Hauck, M., Mohr, H., Peterseim, D.: A simple collocation-type approach to numerical stochastic homogenization. arXiv preprint arXiv:2404.01732 (2024)
40. Hauck M Peterseim D Super-localization of elliptic multiscale problems Math. Comp. 2023 92 341 981 1003
Hauck, M., Peterseim, D.: Super-localization of elliptic multiscale problems. Math. Comp. 92(341), 981–1003 (2023)
41. Hughes TJ Franca LP Hulbert GM A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations Comput. Methods Appl. Mech. Eng. 1989 73 2 173 189
Hughes, T.J., Franca, L.P., Hulbert, G.M.: A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations. Comput. Methods Appl. Mech. Eng. 73(2), 173–189 (1989)
42. Hughes TJR Feijóo GR Mazzei L Quincy JB The variational multiscale method–a paradigm for computational mechanics Comput. Methods Appl. Mech. Engrg. 1998 166 1–2 3 24
Hughes, T.J.R., Feijóo, G.R., Mazzei, L., Quincy, J.B.: The variational multiscale method–a paradigm for computational mechanics. Comput. Methods Appl. Mech. Engrg. 166(1–2), 3–24 (1998)
43. Hughes TJR Sangalli G Variational multiscale analysis: the fine-scale Green’s function, projection, optimization, localization, and stabilized methods SIAM J. Numer. Anal. 2007 45 2 539 557
Hughes, T.J.R., Sangalli, G.: Variational multiscale analysis: the fine-scale Green’s function, projection, optimization, localization, and stabilized methods. SIAM J. Numer. Anal. 45(2), 539–557 (2007)
44. John V Kaya S Layton W A two-level variational multiscale method for convection-dominated convection-diffusion equations Comput. Methods Appl. Mech. Eng. 2006 195 33 4594 4603
John, V., Kaya, S., Layton, W.: A two-level variational multiscale method for convection-dominated convection-diffusion equations. Comput. Methods Appl. Mech. Eng. 195(33), 4594–4603 (2006)
45. Johnson, C.: Numerical solution of partial differential equations by the finite element method. Dover Publications, Inc., Mineola, NY (2009). Reprint of the 1987 edition
46. Kim MY Wheeler MF A multiscale discontinuous Galerkin method for convection-diffusion-reaction problems Comput. Math. Appl. 2014 68 12, Part B 2251 2261
Kim, M.Y., Wheeler, M.F.: A multiscale discontinuous Galerkin method for convection-diffusion-reaction problems. Comput. Math. Appl. 68(12, Part B), 2251–2261 (2014)
47. Kornhuber R Peterseim D Yserentant H An analysis of a class of variational multiscale methods based on subspace decomposition Math. Comp. 2018 87 314 2765 2774
Kornhuber, R., Peterseim, D., Yserentant, H.: An analysis of a class of variational multiscale methods based on subspace decomposition. Math. Comp. 87(314), 2765–2774 (2018)
48. Kuzmin D Algebraic flux correction for finite element discretizations of coupled systems 2007 CIMNE, Barcelona Computational Methods for Coupled Problems in Science and Engineering II 653 656
Kuzmin, D.: Algebraic flux correction for finite element discretizations of coupled systems, pp. 653–656. Computational Methods for Coupled Problems in Science and Engineering II, CIMNE, Barcelona (2007)
49. Larson MG Målqvist A An adaptive variational multiscale method for convection-diffusion problems Commun. Numer. Methods Eng. 2009 25 1 65 79
Larson, M.G., Målqvist, A.: An adaptive variational multiscale method for convection-diffusion problems. Commun. Numer. Methods Eng. 25(1), 65–79 (2009)
50. Li G Peterseim D Schedensack M Error analysis of a variational multiscale stabilization for convection-dominated diffusion equations in two dimensions IMA J. Numer. Anal. 2017 38 3 1229 1253
Li, G., Peterseim, D., Schedensack, M.: Error analysis of a variational multiscale stabilization for convection-dominated diffusion equations in two dimensions. IMA J. Numer. Anal. 38(3), 1229–1253 (2017)
51. Li J Demkowicz L An Lp-DPG method for the convection-diffusion problem Comput. Math. Appl. 2021 95 172 185
Li, J., Demkowicz, L.: An Lp-DPG method for the convection-diffusion problem. Comput. Math. Appl. 95, 172–185 (2021)
52. Målqvist A Peterseim D Localization of elliptic multiscale problems Math. Comp. 2014 83 290 2583 2603
Målqvist, A., Peterseim, D.: Localization of elliptic multiscale problems. Math. Comp. 83(290), 2583–2603 (2014)
53. Målqvist A Peterseim D Numerical homogenization by localized orthogonal decomposition, SIAM Spotlights 2021 Philadelphia, PA Society for Industrial and Applied Mathematics (SIAM)
Målqvist, A., Peterseim, D.: Numerical homogenization by localized orthogonal decomposition, SIAM Spotlights, vol. 5. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA (2021)
54. Melenk, J.M.: Hp-finite element methods for singular perturbations. Lecture Notes Math. 1796 (2002)
55. Miller, J.J.H., O’Riordan, E., Shishkin, G.I.: Fitted numerical methods for singular perturbation problems, revised edn. World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ (2012). 10.1142/9789814390743. Error estimates in the maximum norm for linear problems in one and two dimensions
56. Målqvist A Multiscale methods for elliptic problems Multiscale Model. Simul. 2011 9 3 1064 1086
Målqvist, A.: Multiscale methods for elliptic problems. Multiscale Model. Simul. 9(3), 1064–1086 (2011)
57. Owhadi H Scovel C Operator-adapted wavelets, fast solvers, and numerical homogenization, Cambridge Monographs on Applied and Computational Mathematics 2019 Cambridge Cambridge University Press
Owhadi, H., Scovel, C.: Operator-adapted wavelets, fast solvers, and numerical homogenization, Cambridge Monographs on Applied and Computational Mathematics, vol. 35. Cambridge University Press, Cambridge (2019)
58. Park PJ Hou TY Multiscale numerical methods for singularly perturbed convection-diffusion equations Int. J. Comput. Methods 2004 01 01 17 65
Park, P.J., Hou, T.Y.: Multiscale numerical methods for singularly perturbed convection-diffusion equations. Int. J. Comput. Methods 01(01), 17–65 (2004)
59. Payne LE Weinberger HF An optimal Poincaré inequality for convex domains Arch. Ration. Mech. Anal. 1960 5 1 286 292
Payne, L.E., Weinberger, H.F.: An optimal Poincaré inequality for convex domains. Arch. Ration. Mech. Anal. 5(1), 286–292 (1960)
60. Peterseim, D.: Variational multiscale stabilization and the exponential decay of fine-scale correctors. In: Building bridges: connections and challenges in modern approaches to numerical partial differential equations, Lect. Notes Comput. Sci. Eng., vol. 114, pp. 341–367. Springer, Cham (2016)
61. Qiu W Shi K An HDG method for convection diffusion equation J. Sci. Comput. 2016 66 1 346 357
Qiu, W., Shi, K.: An HDG method for convection diffusion equation. J. Sci. Comput. 66(1), 346–357 (2016)
62. Roos HG Stynes M Tobiska L Robust numerical methods for singularly perturbed differential equations, Springer Series in Computational Mathematics 2008 2 Berlin, Berlin Springer-Verlag
Roos, H.G., Stynes, M., Tobiska, L.: Robust numerical methods for singularly perturbed differential equations, Springer Series in Computational Mathematics, vol. 24, 2nd edn. Springer-Verlag, Berlin, Berlin (2008)
63. Simon K Behrens J Semi-lagrangian subgrid reconstruction for advection-dominant multiscale problems with rough data J. Sci. Comput. 2021 87 2 1 33
Simon, K., Behrens, J.: Semi-lagrangian subgrid reconstruction for advection-dominant multiscale problems with rough data. J. Sci. Comput. 87(2), 1–33 (2021)
64. Xie C Wang G Feng X Variational multiscale virtual element method for the convection-dominated diffusion problem Appl. Math. Lett. 2021 117 107077
Xie, C., Wang, G., Feng, X.: Variational multiscale virtual element method for the convection-dominated diffusion problem. Appl. Math. Lett. 117, 107077 (2021)
65. Zhao L Chung E Constraint energy minimizing generalized multiscale finite element method for convection diffusion equation Multiscale Mod. Simul. 2023 21 2 735 752
Zhao, L., Chung, E.: Constraint energy minimizing generalized multiscale finite element method for convection diffusion equation. Multiscale Mod. Simul. 21(2), 735–752 (2023)
