{'fraction_non_alphanumeric': 0.04413357585170955, 'fraction_numerical': 0.014101489853058388, 'mean_word_length': 4.272222820816722, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Localized model reduction for parameterized problems

In this contribution we present a survey of concepts in localized model order reduction methods for parameterized partial differential equations. The key concept of localized model order reduction is to construct local reduced spaces that have only support on part of the domain and compute a global approximation by a suitable coupling of the local spaces. In detail, we show how optimal local approximation spaces can be constructed and approximated by random sampling. An overview of possible conforming and non-conforming couplings of the local spaces is provided and corresponding localized a posteriori error estimates are derived. We introduce concepts of local basis enrichment, which includes a discussion of adaptivity. Implementational aspects of localized model reduction methods are addressed. Finally, we illustrate the presented concepts for multiscale, linear elasticity and fluid-flow problems, providing several numerical experiments. This work has been accepted as a chapter in P. Benner, S. Grivet-Talocia, A. Quarteroni, G. Rozza, W.H.A. Schilders, L.M. Sileira. Handbook on Model Order Reduction. Walter De Gruyter GmbH, Berlin, 2019+.


Introduction
Projection based model order reduction has become a mature technique for simulations of large classes of parameterized systems; for an introduction, we refer to the text books and survey [BCOW17,HRS16,QN16,BGW15]. However, especially for large-scale and multi-scale problems the "standard" model order reduction approach exhibits several limitations: Curse of parameter dimensionality in the sense that many parameters require prohibitively large reduced spaces, no topological flexibility, and possibly high computational costs and storage requirements in the offline stage for instance due to large computational domains. Localized model order reduction methods, which combine approaches from model order reduction, multiscale methods and/or domain decomposition techniques, overcome or significantly mitigate those limitations. As an further advantage, they allow using reduced spaces of different dimensions in different parts of the computational domain and accommodate (local) changes of the geometry and the partial differential equation (PDE) in the online stage. The key idea of localized model order reduction is to construct local reduced spaces on (unions of) subdomains of the decomposed computational domain and couple the local reduced spaces across interfaces either in a conforming or nonconforming manner. In this chapter we investigate localized model order reduction for linear coercive elliptic parameterized problems; inf-sup stable problems have for instance been considered in [HKP13a] and parabolic and nonlinear problems will be briefly dicussed at the end of this chapter.
A combination of domain decomposition and reduced basis (RB) methods has first been considered in the reduced basis element method (RBEM) [MR02,MR04,LMR06], where the local RB approximations are coupled by Lagrange multipliers in a non-conforming manner. The reduced basis hybrid method [IQR12] extends the RBEM by additionally considering a coarse FE discretization on the whole domain to account for continuity of normal stresses in the context of Stokes equations. Alternatively, a non-conforming coupling can be realized say by penalization as in the local reduced basis discontinuous Galerkin approach [KOH11], the localized reduced basis multiscale method (LRBMS) [AHKO12,OS15,ORS17], the discontinuous Galerkin reduced basis element method [APQ16], or the generalized multiscale discontinuous Galerkin method [CEL17a]. The static condensation reduced basis element (scRBE) method [HKP13b,HKP13a,EP13,SP16] combines intra-element RB approximations similar to the RBEM with coupling techniques from CMS resulting in a conforming approximation. A similar approach is pursued by the ArbiLoMod [BEOR17] that also allows for arbitrary (non-parametric) local changes of the underlying equations and/or the geometry.
In the context of the proper generalized decomposition (PGD) method (for a review see for instance [CA10,CLC11,CKL14]) a domain decomoposition strategy has been proposed in [HNC18] and in [PEV10] hierarchical model reduction [VB81,PEV10,OS14b,SO17] has been combined with an iterative substructuring method.
Concerning the generation of local approximation spaces we focus on empirical training (see for instance [EP13,BL11,SP16]), i.e. local reduced spaces generated from local solutions of the PDE, and adaptive basis enrichment. In detail, we present local approximation spaces that are optimal in the sense of Kolmogorov and can be constructed by solving a local so-called transfer eigenvalue problem on the space of local solutions of the PDE. Optimal local approximation spaces for subdomains have first been proposed in [BL11] and for interfaces and parametrized PDEs in [SP16]. We will also show how those optimal approximation spaces can approximated by random sampling [BS18].
A localizable a posteriori error estimator is crucial for an adaptive enrichment of the local reduced spaces where the reduced approximation is not accurate enough. Such an adaptive basis enrichment is one way to approach "optimal" computational complexity within outer-loop applications such as optimal control, inverse problems or Monte Carlo methods. With this respect, we will also present a framework for localized residual based error control [BEOR17,Sme15] as well as localized a posteriori error estimation based on flux reconstruction [ESV10,OS15].
Naturally, the presented methods for localized model reduction share a lot of features with domain decomposition techniques and multiscale methods. We particularly refer to domain decomposition and preconditioning techniques with multiscale coarse spaces such as [AH02,GLS07,GE10] or the more recent contributions [SDH + 14,GL17,HKKR18]. In the context of the FETI-DP iterative substructuring method we refer to [MS07,KRR16]. For multiscale problems there has been a tremendous development of suitable numerical methods in the last two decades including the multiscale finite element method (MsFEM) [HW97,EHG04,EH09,HOS14], the heterogeneous multiscale method (HMM) [EE03,EE05,Ohl05,Abd05], the variational multiscale method (VMM) [Hug95,HFMQ98,LM05] or the more recent local orthogonal decomposition (LOD) [MP14,HMP14]. Model reduction can be used to accelerate the solution of localized problems which occur in multiscale methods, see e.g. [AB12,Abd15]. Similar to the methods presented in this chapter the generalized multiscale finite element method (GMsFEM) [EGH13,CEL14b,CEL18a] relies on a Galerkin projection on local subspaces, but in contrast uses ideas from multiscale methods to construct the local bases. A connection between multiscale methods and domain decomposition has recently been investigated in [KPY18,KPY17,KY16].
This chapter is organized as follows. In Sec. 2 we introduce the problem setting and basic notation for localized model order reduction of coercive variational problems. Concepts for conforming and non-conforming coupling of approximation spaces are presented in Sec. 3. Sec. 4 deals with the preparation of local approximation spaces. Particularly, the construction of optimal local approximation spaces and their approximation via random sampling is presented and illustrated with numerical experiments for linear elasticity. In Sec. 5 we present two abstract frameworks for localized a posteriori error estimation and give exemplifications for conforming and non-conforming localized model reduction approaches. Localized a posteriori error estimators are the key ingredient for basis enrichment strategies and online adaptivity that are presented in Sec. 6. Computational aspects are discussed in Sec. 7 and numerical experiments for multiscale problems and fluid flow are presented in Sec. 8. We conclude by showing possible extensions to parabolic and nonlinear problems in Sec. 9.

Parameterized partial differential equations and localization
Let Ω ⊂ R d , d = 1, 2, 3, be a large, bounded domain with Lipschitz boundary. Let us further introduce a Hilbert space V such that [H 1 0 (Ω)] z ⊂ V ⊂ [H 1 (Ω)] z , z = 1, . . . , d and denote by V the dual space of V . Moreover, we introduce the compact set of admissible parameters P ⊂ R p , p ∈ N. We consider the following variational problem.
Definition 2.1 (Parameterized coercive problem in variational form). For any parameter µ ∈ P find u(µ) ∈ V , such that Here, f (·; µ) ∈ V and a(·, ·; µ) : V × V → R denote parametric linear and bilinear forms, the latter being continuous and coercive w.r.t. the norm · V induced by the inner product (·, ·) V : V × V → R. That is, there exist constants 0 < α ≤ α(µ) ≤ γ(µ) ≤ γ, such that for any µ ∈ P, Let us denote the energy norm of u for parameter µ as |||u||| µ := a(u, u; µ) 1/2 . Problem (1) thus admits a unique solution for all µ ∈ P owing to the Lax-Milgram theorem. Examples for (1) include elliptic multiscale problems, incompressible fluid flow or linear elasticity as detailed in the following. We will consider Neumann boundary conditions on Γ N and Dirichlet boundary conditions on Γ D , where Γ N , Γ D are non-overlapping and Γ N ∪ Γ D = ∂Ω. To simplify notations, homogenous boundary conditions on ∂Ω will be prescribed in most places.
Example 2.2 (Parametric elliptic multiscale problems). With V = H 1 0 (Ω), the pressure equation in the context of two-phase flow in porous media (obtained through Darcy's law) reads: given a collection of sources and sinks q ∈ L 2 (Ω), a parametric and possibly highly heterogeneous permeability field κ : P → L ∞ (Ω) d×d , find for each µ ∈ P the global pressure u(µ) ∈ V , such that −∇ · κ(µ)∇u(µ) = q in a weak sense in V .

(2)
If the smallest eigenvalue of κ(µ) is bounded from below away from zero for all µ ∈ P, we can consider this to be an example of Definition 2.1 by setting a(u, v; µ) := Ω (κ(µ)∇u) · ∇v dx and f (v; µ) := Ω qv dx. In the context of instationary twophase flow, (2) needs to be solved in each time step for varying total mobilities (modelled by the parametric nature of κ), while the permeability field κ typically resolves fine geological structures and thus requires a very fine computational grid compared to the size of Ω (see [OS15] and the references therein). In the steady case it can be formulated as follows: (3) where (y, p) are the velocity and the pressure fields defined on the computational domain Ω. The first equation expresses the linear momentum conservation, the second one the mass conservation, which is also called the continuity equation. Here f denotes a forcing term per unit mass, g D and g N are the functions addressing the Dirichlet and Neumann boundary conditions respectively on Γ D and Γ N . The parameter ν = σ/ρ denotes the kinematic viscosity, being ρ the density and σ the viscosity of the fluid. Navier-Stokes equations correspond to the case δ = 1, here we consider only δ = 0, the convective term is neglected, obtaining the steady Stokes equations, which provide a model in the case of slow motion of fluids with very high viscosity.
We denote the functional spaces for velocity and pressure fields by X = (H 1 0,Γ D (Ω)) d , Q = L 2 (Ω), respectively, where H 1 0,Γ D (Ω) = {y ∈ H 1 (Ω) : y| Γ D = 0}. Moreover, for simplicity, we assume that g D = 0 (otherwise the lift function is required). The corresponding weak form of the Stokes equations (3) reads: find (y, p) ∈ X × Q such that In a parameterized setting, the input-parameter vector µ may characterize either the geometrical configuration or physical properties, boundary data and sources of the problem. Denoting by V the product space given by V = X × Q, defining by u(µ) = (y(µ), p(µ)) ∈ V and v = (w, q), the parametrized abstract formulation (2.3) can be rewritten in the following form: find u(µ) = (y(µ), p(µ)) ∈ V s.t.
To obtain approximate solutions of (1) we presume we have an appropriate grid-based numerical method at hand (the full order model, FOM), yielding a high (but finite) dimensional approximation space V h . We consider conforming continuous Galerkin finite elements (FE), where V h ⊂ V , and nonconforming discontinuous Galerkin (DG) or finite volume (FV) schemes, where V h ⊂ V (in which case we require broken Sobolev spaces for our analysis, see Section 3.2.2). As a starting point for localized model reduction we require the FOM space to be decomposable into "local" spaces, which we will make more precise shortly. While a localizing space decomposition could in general stem from any clustering of the degrees of freedom (DoF) of V h (see for instance [Car15]), we are particularly interested in local approximation spaces which are associated with a domain decomposition of the physical domain.  We can thus think of the domain decomposition as a usual grid, but without the requirements of T H to actually resolve any data functions of the PDE. Given such a domain decomposition, we can abstractly define a localizing space decomposition.
Definition 2.6 (Localizing space decomposition). Let the FOM space V h be a finite dimensional Hilbert space with inner product and induced norm · 2 V h = (·, ·) V h . We call the direct sum decomposition of V h into subdomain spaces, interface spaces, edge spaces and vertex spaces, a localizing space decomposition.
Note that such a decomposition is not unique and can always be found. Since the reduced space shall inherit this localizing decomposition, its purpose will be threefold: (i) offline, it allows for an independent and localized generation of the local reduced approximation spaces (compare Section 4), (ii) it allows to define and alter the physical domain Ω online, given that local approximation spaces for certain reference subdomains have been prepared offline, and (iii) it allows to adapt a local approximation space online (by adding basis functions or changing the local grid), while only requiring an update of local and neighboring prepared quantities (compare Section 6). For actual examples of space decompositions we refer to Section 3.
Abstractly, we do not impose any further assumptions on the FOM as well as the reduced order model (ROM). However, given the (bi-)linearity of a and f , the computational benefits of the localizing space decomposition are apparent (and are made more precise throughout the rest of this chapter). Since we allow for non conforming approximations, in general we need to consider discrete counterparts of a and f which are only defined on the FOM space V h and not necessarily on V , where we again refer to the following sections for examples.
Definition 2.7 (Locally decomposed full order model (FOM)). Let V h be locally decomposable as in definition (9), and let a h (·, ·; µ) : V h × V h → R and f h (·; µ) ∈ V h denote discrete variants of a and f , respectively, which are continuous and coercive w.r.t. the inner product of V h . For each µ ∈ P, find u h (µ) ∈ V h such that The idea of projection-based localized model order reduction is to consider a local reduced approximation space for each element of the localizing space decomposition (9), in order to obtain a similarly decomposed reduced space V N ⊂ V h : Similar to standard projection based model order reduction, we obtain the ROM simply by Galerkin projection of the locally decomposed FOM (10) onto this locally decomposed reduced space.
Definition 2.8 (Locally decomposed reduced order model (ROM)). Given a locally decomposed reduced space as in (11), for each µ ∈ P, find u N (µ) ∈ V N such that The main questions remain: (i) how to choose good local reduced approximation spaces to guarantee accurate and at the same time efficient reduced order approximations, (ii) how to benefit from the localization of V N , that is how to obtain an offline/online decomposed scheme and in particular how to couple these local reduced approximation spaces, and (iii) how to adaptively enrich these local reduced approximation spaces online, if required. These topics will be answered throughout the remainder of this chapter, starting with examples of how to obtain localized FOMs from standard discretization schemes and how to couple the resulting local reduced spaces. Therefore, we introduce local grids τ h (Ω m ) on each subdomain Ω m ⊂ T H , which we presume to resolve all data functions of the underlying PDE. As an analytical tool, we also define the global fine grid by τ h = ∪ Ωm∈T H τ h (Ω m ), which is usually not required in practical computations. For simplicity, we require the local grids of two subdomains Ω m , Ω m ∈ T H to match along the shared interface γ m,m ∈T γ H and denote by τ γ h (γ m,m ) the corresponding set of all facets of τ h which lie on γ m,m . Finally, we require that Γ does not cut any grid cells.
3. Coupling local approximation spaces 3.1. Conforming approach. There are various ways to couple local reduced spaces such that we obtain a conforming approximation, such as the partition of unity method [BM97] or the generalized finite element method (GFEM) [BCO94,BBO04,BM97,BL11]. However, in this section we focus on the decomposition into interface spaces and intra-element spaces, where the coupling is performed via the coupling or interface modes spanning the interface space.
The formulation (13) can then be used to derive an equation that solely acts on functions on the interface but nevertheless uniquely defines the solution u(µ) of (1). To that end, we introduce a parameter-dependent lifting operator E Γ→Ω (µ) : Then, we can rewrite the solution u(µ) as Inserting (15) into (13c) and choosing E m = E Γ→Ωm (µ) yields the Steklov-Poincaré interface equation: Let us notice that the Steklov-Poincaré interface equation and its discrete, algebraic analogon, the Schur complement system, are at the base of iterative substructuring methods (see [QV05,TW05]), which have been combined with the reduced basis method in [MH14]. We may finally define a space associated with the interface V Γ as V Γ = {E Γ→Ω (µ)ζ ∈ V : ζ ∈ Λ} and obtain the decomposition This decomposition is a-orthogonal thanks to (14).
While the computation of the (harmonic) lifting operators is inherently local (see (14)), the Steklov-Poincaré interface equation is posed on the whole interface Γ. To localize the latter we decompose V Γ as we will describe next.
3.1.2. A conforming, localized reduced order approximation. First, we determine basis functions associated with the vertices v ∈ T v H . One common approach [HL10,JBL11,BEOR17] is to require that a basis function Here, x v are the (global) coordinates of the vertex v and V m h;0 := {v ∈ V m h : v = 0 on ∂Ω m \ Γ N } . To uniquely define ψ v we need to prescribe the respective values on Γ. We may for instance require ψ v to be linear on the respective edges or bilinear on the respective interfaces (see e. g. [BEOR17]). For multiscale problems in two space dimensions with a permeability κ(x 1 , x 2 ;μ) it has been suggested in [HW97] to prescribe

}.
Similarly as for the vertices we may then define associated basis functions that have support on the union of subdomains that share the respective edge or interface: such that (21) (ψ e k , w) V m = 0 ∀w ∈ V m h;0 and ψ e k | e = χ e k . Again, we need to provide the value of ψ e k on the interfaces sharing the edge e ∈ T e H in order to uniquely define ψ e k . Similarly to above we may require that ψ e k is linear on the respective interfaces as suggested for instance in [BEOR17] or define a function which takes into account also the coefficient function.
Note that if the interfaces are mutually disjoint, which is for instance the case if we associate the subdomains Ω m , m = 1, . . For N γ N γ h;0 and N e N e h;0 may now define the reduced space associated with Γ as follows: Such reduced interface spaces are for instance employed in (adaptive) component mode synthesis (CMS) [HL10,JBL11], the static condensation reduced basis element (scRBE) method [HKP13b, HKP13a, EP13, SP16] for mutually disjoint interfaces, or in the ArbiLoMod [BEOR17]. Subspaces of V N Γ are considered in certain multiscale methods. For example in the MsFEM of Hou and Wu [HW97] the reduced space is spanned by the basis functions ψ v , v ∈ T v H . For further relations between CMS, MsFEM and GFEM we refer e.g. to [HL10].
Recall that the basis functions associated with the vertices, edges, and interfaces have all been computed w.r.t. an inner product that does not depend on the parameter (see (18), (20), (21)). Therefore, we finally assume that we have also given reduced spaces V m N ;0 := span{ζ m 1 , . . . , ζ m N m } ⊂ V m h;0 , m = 1, . . . , M , that will account for parameter variations. In detail, we obtain approximationsψ * k (µ), * = v, e, γ by solving Note, that bothb * k (µ), * = v, e, γ andb m (µ) can be interpreted as intra-element RB approximations. The corresponding reduced spaces V m N ;0 , m = 1, . . . , M can for instance be constructed from h;0 , respectively, via a standard greedy algorithm or a POD. 2 Let us also remark that for instance in the scRBE method for the approximation of each basis function ψ * k , * = v, e, γ a different RB space is considered, to further reduce the size of problems (23), (24). Finally, we define the reduced spaces The global reduced approximation u N (µ) can then be computed by performing a Galerkin projection onto the reduced space V Γ N (µ) or a Petrov-Galerkin approximation using V Γ N (µ) as a trial and V Γ N as a test space (see e.g. [EP14,Sme15]). Instead of eliminating the volume degrees of freedom via (23), (24), u N (µ) can also directly be determined by performing a Galerkin projection onto V N (see for instance [BEOR17]). Similarly, for CMS and a fixed parameter a Galerkin projection onto V N may be performed to compute the reduced solution; here, the reduced space V m N ;0 is constructed from an eigenvalue problem and does not account for parameter variations (see e.g. [HL10]). Finally, in the reduced basis-domain decomposition-finite element (RDF) method [IQR16] the reduced space V N is chosen as a direct sum of M m=1 V m N ;0 and standard FE spaces defined on the interface or on a (small) area around the interface. Here, the intra-element reduced spaces V m N ;0 are constructed via a greedy algorithm considering a parametrized linear combination of standard Lagrange basis functions or Fourier modes as boundary conditions. Then, a Galerkin projection on V N is performed to compute u N (µ).
3.2. Non-conforming approach. With the term non-conforming approach we want to classify a set of alternative techniques to solve the reduced problem on the global computational domain. A first approach consists in considering a global system of equations given by local parametrized problems and additional equations ensuring the matching between the different subdomains through the use of Lagrange multipliers. This approach has been used for solving elliptic equation in [MR02,MR04] and Stokes equations in [LMR06,IQR12].
Another approach consists in coupling local FOM spaces by interior penalty bilinear forms, inspired by discontinuous Galerkin FEM. Here, we refer to the discontinuous Galerkin reduced basis element method [CHM11, APQ16, PGQ16] and the local reduced basis discontinuous Galerkin approach [KOH11]. A discontinuous Galerkin approach with local POD modes has been presented in [FIL18]. In the context of multiscale problems (cf. Example 2.2), the generalized multiscale discontinuous Galerkin method has been proposed in [CEL17a,CEL18b] and used 2 Note that in actual practice one would construct the reduced bases only on a certain number < M of reference domains; see for instance [HKP13b].
for solving the heat problem with phase change in [SVV18]. In this chapter we are going to present the localized reduced basis method (LRBMS) in Subsection 3.2.2 below. LRBMS has been introduced in [AHKO12] and analyzed in [OS14a,OS15] for elliptic and in [ORS17] for parabolic problems. Applications to the simulation of two phase flow in porous media have been addressed in [KFH + 15] and to battery simulation with resolved electrodes in [OR17].
3.2.1. Non-conforming coupling based on Lagrange multipliers. We want to reformulate the problem (13), with the idea that exact coincidence of the traces of the discrete functions (equation (13b)) is generally too stringent, and may, in fact, lead to imposing u m = 0 on the internal interfaces; thus, the gluing process can be done in a dual way through Lagrange multipliers. We assume that local basis functions are computed in each subdomain Ω m , m = 1, . . . , M by solving local parametrized variational problems coming from the original problem (1) with proper boundary conditions along the boundaries which correspond to internal ones in the original domain. The choice of the boundary conditions is strongly related to the problem aimed to be solved. Thus, local reduced spaces are defined via these local solutions and denoted by V m N , m = 1, . . . , M . Possible ways to construct V m N are presented in section 4. If two or more subdomains are characterized by the same type of parameter and the same type of boundary conditions, the same local reduced space can be associated to those subdomains. For simplicity we consider different spaces for each different subdomain. We define the following operator: where m, m ∈ {1, · · · , M }, Γ m,m is the interface between two adjacent sub-domains denoted with the indices m and m respectively, and W m,m is the Lagrange multiplier space defined on this interface. Typical choices for the latter are low-order polynomial spaces [MR02,IQR12] or spaces constructed from snapshots (and their derivatives) [MR04].
A basis for W m,m can then for instance be provided by the characteristic Lagrange polynomials ψ q , q = 1, · · · , Q m,m associated with the Q m,m nodes of Γ m,m . If we suppose that Ω has M − 1 internal interfaces, Γ m,m+1 , m = 1, · · · , M − 1, the reduced global problem of this approach reads:

3.2.2.
Non-conforming coupling based on interior penalties. We demonstrate how to obtain a localized FOM by applying ideas from interior penalty (IP) DG schemes w.r.t. the domain decomposition T H in the context of the parametric multiscale Example 2.2. To define the localized FOM, we presume we are given a discretization on the full global grid τ h (which is not used in actual computations), which we make precise by specifying the approximation space with an associated inner product and discrete variants of a and f . As a common ground for the analysis of conforming as well as non conforming schemes, we introduce the broken Sobolev space , for a given grid τ h (ω) of some domain ω ⊆ Ω and s ≥ 1, and associated broken gradient operator ∇ h : Example 3.1 (Continuous Galerkin (CG) FEM). The CG FEM scheme for the conforming approximation of Example 2.2 w.r.t. the full global grid τ h is given by the conforming approximation space of order k ≥ 1, (again noting that their respective restriction to V coincide with a and f ).
The definition of the non conforming scheme is more involved. We denote the set of faces of τ h by τ γ h and to each face σ ∈ τ γ h , we assign a unique normal n σ ∈ R d pointing away from t + , where the face may be either an inner face σ ∈ τ γ h , given by the intersection of two grid elements t + , t − ∈ τ h , σ = t + ∩ t − , or a boundary face σ ∈ τ γ h , given by σ = t + ∩ ∂Ω for some t + ∈ τ h . Since functions in the broken Sobolev space are two-valued on grid faces, we introduce the mean · and jump [·] on a boundary face by v := [v] := v| t + and by v := 1 2 (v| t + + v| t − ) and [v] := v| t + − v| t − , respectively, on any other face.
Considering the family of interior penalty discontinuous Galerkin (DG) schemes, we present the symmetric variant for ease of notation, and refer to the symmetric weighted variant [ESZ09], which is particularly well suited for multi-scale problems with highly varying or high-contrast coefficients.
Example 3.2 (Interior penalty (IP) discontinuous Galerkin (DG) FEM). The symmetric IPDG FEM scheme for the non conforming approximation of Example 2.2 w.r.t. the full global grid τ h is given by the non conforming approximation space of order k ≥ 1, with the face bilinear form a σ for any σ ∈ τ γ h given by ds and a user-dependent penalty weight w σ > 0, such that a DG h is continuous and coercive w.r.t. the above inner product.
The main idea of an IP localized FOM is to consider the restriction of either of the above discretization schemes to each subdomain of the domain decomposition, and to again couple those with IP techniques along the interfaces of the subdomain. We thus choose * ∈ {CG, DG} and obtain the localized FOM space in the sense of Definition 2.6 as a direct sum of subdomain spaces (with empty interface, edge and vertex spaces) We also define the linear functional f h : V h → R by f h := f * h and, in a similar manner as above, the non conforming bilinear form a h (·, ·; µ) : We have thus fully specified a localized FOM in the sense of Definition 2.7 and comment on two special cases: for * = CG and a trivial domain decomposition of a single subdomain, T H = {Ω}, we obtain the above standard CG FEM while for * = DG the resulting FOM coincides with the above standard symmetric IPDG FEM.
To make the coupling more precise, we may rearrange the above terms to obtain a localization of a h w.r.t. the domain decomposition in the sense of with the subdomain and interface bilinear forms respectively, for all 1 ≤ m ≤ M and all Γ ∈T γ H , with the subdomains Ω + , Ω − ∈ T H sharing the interface Γ . Now given a local reduced space V m N ⊂ V m h for each subdomain we obtain the locally decomposed broken reduced space in the sense of (11) by Using the above decomposition of a h into subdomain and interface contributions, we can readily observe that the locally decomposed ROM can be offline/online decomposed by local computations: namely by projection of the subdomain bilinear forms a m h (·, ·; µ) onto V m N × V m N and the interface bilinear forms a Γ h (·, ·; µ) onto V m N × V n N , with 1 ≤ m, n ≤ M , such that Ω + = Ω m and Ω − = Ω n , respectively. We thus obtain a sparse matrix representation of the resulting reduced system, with a sparsity pattern which coincides with the one from standard IPDG schemes.

Preparation of local approximation spaces
Both, couplings that yield a conforming and non-conforming approximation require either reduced spaces Λ γ N γ ⊂ V h | γ for interfaces and/or edges Λ e N e ⊂ V h | e (see subsection 3.1) or reduced spaces V m N (see subsection 3.2) or both. As the generation of edge basis functions can be done analogously to the construction of interface basis functions we restrict ourselves to the latter in order to simplify notation. To fix the setting we thus consider the task of finding a suitable reduced space either on a subdomain Ω m Possible geometric configurations of the oversampling domain Ω out are illustrated in Fig. 1.
We will first briefly discuss in subsection 4.1 reduced spaces that are spanned by polynomials or solutions of "standard" eigenvalue problems and are thus related to the spectral element method or hp-FEM. Subsequently, in subsection 4.2 we will present reduced spaces that are generated from local solutions of the PDE, are thus of empirical nature, and are optimal in the sense of Kolmogorov. We will also show how those optimal basis functions can be efficiently and accurately approximated by means of random sampling. A combination of domain decomposition and RB methods has first been considered in the reduced basis element method (RBEM) [MR02]. Here, inspired by the the mortar spectral element method [BMP94], the Lagrange multiplier space W m,m as defined in Subsection 3.2 is chosen as a low-order polynomial space. The reduced basis hybrid method [IQR12] extends the RBEM by additionally considering a coarse FE discretization on the whole domain to account for continuity of normal stresses and also employs a low-order polynomial Lagrange multiplier space on the interface. For the scRBE method a reduced interface space spanned by the eigenvectors of a discrete generalized eigenvalue problem based on the Laplacian has been suggested in [HKP13b,HKP13a] and eigenmodes of a singular Sturm-Liouville eigenproblem have been used in [EP13]. Finally, in the RDF method [IQR16] a standard FE space is considered on the interface or on a (small) area around the interface.

4.2.
Local approximation spaces based on empirical training. In this subsection we are concerned with local approximation spaces that are constructed from local solutions of the PDE; those approaches are often called empirical. Basis functions on the interfaces selected from local snapshots are for instance suggested in [EP13], where an empirical pairwise training procedure for interface reduction within the scRBE context is developed, and within a heterogeneous domain decomposition method in [MRH15]. Local approximation spaces that are optimal in the sense of Kolmogorov have been introduced for subdomains within the GFEM in [BL11] for parameter-independent PDEs and for interfaces within static condensation procedures [SP16] for parametrized PDEs. While the authors of [SP16] introduce and analyze a spectral greedy algorithm to deal with parameter variations, [TP18] suggests using a POD making use of the hierarchical approximate POD [HLR18]. Those optimal local spaces both allow for a rigorous a priori theory and yield a rapidly (and often exponentially) convergent approximation; in certain cases the superalgebraic convergence can be proved [BL11]. Recently, in [TP18,Tad16] the results in [BL11, BHL14, SP16] have been generalized from linear differential operators whose associated bilinear form is coercive to elliptic, inf-sup stable ones. In [BS18] it has been shown that those optimal local approximation spaces can be efficiently approximated by transferring methods from randomized numerical linear algebra [HMT11]; the local approximation approximation spaces in [BS18] are constructed from local solutions of the PDE with random boundary conditions. Local reduced spaces generated from random snapshots have also been suggested in [BEOR17,EP13] and methods from randomized linear algebra have been exploited in the FETI-2λ domain decomposition method in [WV15] and in [CEGL16] for the generalized multiscale finite element method.
We will first present the optimal local approximation spaces as introduced in [BL11, SP16] for a fixed parameter P = {μ}, subsequently discuss their approximation via random sampling, and conclude this subsection with the discussion of the general case P = {μ}. To simplify notation we will omitμ as long as it is fixed. 4.2.1. Optimal local approximation spaces for P = {μ}. To enable maximum flexibility regarding the shape of Ω on the user's side, we assume that we do not have any a priori knowledge of the shape of Ω when constructing the reduced order model. We thus know that the global solution u satisfies the considered PDE locally on Ω out but suppose that the trace of u on ∂Ω out is unknown to us. Therefore, we aim at approximating all local solutions u loc of with arbitrary Dirichlet boundary conditions on Γ out . Here, the Hilbert space V loc is defined such that [H 1 0 (Ω out )] z ⊂ V loc ⊂ [H 1 (Ω out )] z , z = 1, . . . , d, respecting the boundary conditions on ∂Ω, and a loc : R are local bilinear and linear forms. We will first restrict ourselves to the case f loc = 0, g D = 0, and ∂Ω m ∩ Γ D = ∅; the general case will be dealt with at the end of this subsubsection. We may then define the space of all local solutions of the PDE as As suggested in [BL11,SP16] we introduce a transfer operator T : S → R for Hilbert spaces S and R, where S = {w| Γout : w ∈ H}. We define T for interfaces or subdomains, respectively, for w ∈ H as Here, P D , D ⊂ Ω out , denotes an orthogonal projection onto the kernel of the bilinear form; for further details see [BS18,SP16]. In the case of heat conduction we would for instance subtract the mean value of the respective function on D. Note that subtracting this projection is necessary to prove compactness of the transfer operator T . The key argument to show compactness of T is Caccioppoli's inequality, which estimates the energy norm of a function in H on Ω m in terms of the L 2 -norm on Ω out of the respective function. Using the Hilbert-Schmidt theorem and Theorem 2.2 in [Pin85, Chapter 4] it can then be shown that certain eigenfunctions of T * T span the optimal local approximation space, where T * : R → S denotes the adjoint operator of T . As we aim at approximating H and thus a whole set of functions, the concept of optimality of Kolmogorov [Kol36] is used: A subspace R n ⊂ R of dimension at most n for which holds is called an optimal subspace for d n (T (S); R), where the Kolmogorov n-width d n (T (S); R) is defined as We summarize the findings about the optimal local approximation spaces in the following theorem.
Theorem 4.1 (Optimal local approximation spaces [BL11,SP16]). The optimal approximation space for d n (T (S); R) is given by .., n, and λ j are the largest n eigenvalues and φ j the corresponding eigenfunctions that satisfy the transfer eigenvalue problem: Moreover, we have: Remark 4.2. We emphasize that the optimal space R n is optimal in the sense of Kolmogorov for the approximation of the range of T and not necessarily for the approximation of u(µ). Moreover, we remark that χ sp i are the left singular vectors and √ λ i the singular values of T .
Next, for f loc = 0 but still g D = 0 we solve the problem: To take non-homogeneous Dirichlet boundary conditions into account one can proceed for instance with a standard lifting approach, adjusting f loc accordingly. Note that for homogeneous boundary conditions we proceed very similar to above, prescribing "arbitrary" boundary conditions on Γ out and homogeneous boundary conditions on ∂Ω ∩ ∂Ω out . The optimal local approximation space for subdomains are then defined as ) and similarly for interfaces as 0 for the bilinear form a m defined in subsection 3.1. In case ∂Ω m ∩ Γ D = ∅ all modifications in this subsubsection involving the kernel of the bilinear form are waived.
The result in (36) can be exploited to derive an a priori error bound for the approximation error between the solution u(μ) of (1) still for a fixed reference parameterμ and the optimal static condensation approximation u n (μ) as stated in the following proposition: Proposition 4.3 (A priori error bound [SP16]). Assume that the interfaces γ ∈ T γ H are mutually disjoint, that all interfaces have the same geometry, and that each Ω m , m = 1, . . . , M has exactly two interfaces. Let u(μ) be the (exact) solution of (1) for a fixed parameterμ. Moreover, let u n+ (μ) be the static condensation approximation defined in subsection 3.1, where we employ the optimal interface space R + n for each γ ∈ T γ H and assume that the error due to the intra-element RB approximation is zero. Then, we have the following a priori error bound: where #γ denotes the number of interfaces in T γ H and λ γ n+1 is the n + 1-th eigenvalue of (35) for the interface γ ∈ T γ H . The constant C γ depends only on the subdomains that share the interface γ and neither on Ω nor on u(μ).
To define reduced interface spaces Λ γ N γ , γ ∈ T γ H and reduced spaces V m N , m = 1, . . . , M we approximate (35) with finite elements. To that end, we introduce a We may then define the discrete transfer operator T : In order to define a matrix form of the transfer operator we introduce DOF mappings that map the DOFs of S to the DOFs of V h | Ωout and the DOFs of V h | Ωout to the DOFs of R, respectively. Moreover, we introduce the stiffness matrix A loc obtained from the FE discretization of (31), where we assume that in the rows associated with the Dirichlet DOFs the non-diagonal entries are zero and the diagonal entries equal one. By denoting by ζ the FE coefficients of ζ ∈ S and by defining P D as the matrix of the orthogonal projection on the kernel of the bilinear form on D ⊂ Ω out , we obtain the following matrix representation T ∈ R N R ×N S of the transfer operator for subdomains and interfaces Finally, we denote by M S the inner product matrix of S and by M R the inner product matrix of R. Then, the FE approximation of the transfer eigenvalue problem reads as follows: Find the eigenvectors ζ j ∈ R N S and the eigenvalues λ j ∈ R + 0 such that (43) The coefficients of the FE approximation of the basis functions {χ sp h,1 , ..., χ sp h,n } of the discrete optimal local approximation space (44) R n := span{χ sp h,1 , ..., χ sp h,n } are then given by χ sp h,j = T ζ j , j = 1, . . . , n. Adding the representation of the right-hand side, the boundary conditions, and a basis of the kernel of the bilinear form yields the optimal spaces Λ γ N γ and V m N . Note that in actual practice we would not assemble the matrix T. Instead one may solve the PDE locally N S times prescribing the basis functions of S as Dirichlet boundary conditions on Γ out and subsequently assemble and solve the transfer eigenvalue problem. Alternatively, one may pass T implicitly to the Lanczos method. For instance, the implicitly restarted Lanczos method as implemented in ARPACK [LSY98] requires O(n) local solutions of the PDE in each iteration and applications of the adjoint T * . In the next subsubsection we will show how methods from randomized linear algebra [HMT11,DM16,Mah11,MD09] can be used to compute an approximation of the optimal local approximation spaces. However, beforehand, we conclude this subsubsection with some numerical experiments on the transfer eigenvalues and thus via Proposition 4.3 on the convergence behavior of the relative approximation error.
To this end, we present the simplified model for a ship stiffener from [SP16]: We consider Ω out =Ω 1 ∪Ω 2 and Γ m,m = Γ 1,2 =Ω 1 ∩Ω 2 , where Ω 2 is depicted in Fig. 2, Ω 1 is just a shifted version of Ω 2 , and the part of Γ out in Ω 2 is indicated in yellow in Fig. 2. We allow E(µ) to vary in the red areas of the subdomains between 1 and 20 and prescribe E(µ) ≡ 1 in the gray areas; we choose G(µ) = (0, 0, 0) T .
In detail, we consider Ω 1 = (−0.7, 0.   both S and R with a lifting inner product based on the lifting operator E Γ→Ωm (μ) defined in subsubsection 3.1.1; for further details we refer to [SP16]. We consider different values for Young's modulus (ratios) E r i , i = 1, 2 in the red areas of the subdomains and observe in Fig. 3 for the ship stiffener application an exponential convergence of order ≈ e −n of the eigenvalues λ n (µ) and thus the static condensation approximation. We emphasize that we observe in Fig. 3 that the eigenvalues associated with the stiffened plate (E r 1 = E r 2 = 20) decay fastest, while we see the slowest decay for the non-stiffened plate (E r 1 = E r 2 = 1). This is consistent with the expectation that stiffening the plate decreases the deflection of the plate, eliminating the higher eigenmodes. Moreover, an inspection of the optimal interface modes reveals many "classical" mode shapes such as bending or torsional modes of beams and demonstrates again the physical significance of the optimal modes. Also for beams of different shapes, including an I-beam with a crack and thus an irregular domain, an exponential convergence of the transfer eigenvalues and the physical significance of the transfer eigenmodes can be observed; for further details see [SP16]. 4.2.2. Randomized Training. In order to compute an efficient approximation R rand n of R n the adaptive randomized range approximation algorithm 4.1 as suggested in [BS18] iteratively enhances the reduced space with applications of T to a random function until a certain convergence criterion is satisfied.
In detail, in each iteration in line 5 we draw a new random vector r ∈ R N S whose entries are independent and identically distributed random variables with standard normal distribution. Then, we employ the mapping D −1 S : R N S → S to define a unique FE function in S whose coefficients are the components of r.

Algorithm 4.1: Adaptive Randomized Range Approximation
Input : Operator T , target accuracy tol, number of test vectors n t , maximum failure probability ε algofail Output : space R rand n with property P T − P R rand . . , T D −1 S r nt } 2 Compute error estimator factors: Subsequently, we apply the transfer operator T to D −1 S r, meaning that we solve the PDE locally on Ω out with random boundary conditions and restrict the solution to Ω m or Γ m,m ; the resulting function is added to the set of basis functions B. Finally, the basis B is orthonormalized. Note that the orthonormalization is numerically challenging, as the basis functions are nearly linear dependent when span{B} is already a good approximation of the range of T ; in [BS18] using the numerically stable Gram-Schmidt with adaptive re-iteration from [BEOR14] is suggested. The main loop of the algorithm is terminated when the following a posteriori norm estimator is smaller than the desired tolerance tol.
Proposition 4.4 (A probabilistic a posteriori norm estimator [BS18]). Let r i , i = 1, . . . , n t be n t random normal test vectors and λ M S min and λ M S max the smallest and largest eigenvalues of the matrix of the inner product in S. Then, the a posteriori norm estimator (45) ∆(n t , ε testfail ) := c est (n t , ε testfail ) max i∈1,...,nt where c est (n t , ε testfail ) : where the constant c eff (n t , ε testfail ) is defined as and Q −1 is the inverse of the upper normalized incomplete gamma function.
The constant c est (n t , ε testfail ) is calculated in line 3 using N T , which denotes the rank of operator T . In practice N T is unknown and an upper bound for N T such as min(N S , N R ) can be used instead. Note that the term (max t∈M t R ) · c est (n t , ε testfail ) is the norm estimator (45). The test vectors are reused for all iterations.
To finally analyze the failure probability of Algorithm 4.1 we first note that after N T steps we have R rand n = range(T ) and thus T − P R rand n T = 0, yielding the termination of Algorithm 4.1. Using the fact that the a posteriori error estimator defined in (45) is therefore executed at most N T times combined with the probability for one estimate to fail in (46) and an union bound argument we infer that the failure probability for the whole algorithm is ε algofail ≤ N T ε testfail .
Remarkably, the convergence behavior of the reduced space R rand n is only slightly worse than the rate λ n+1 , which is achieved by the optimal local approximation spaces defined in Theorem 4.1: . It can be observed in numerical experiments that the a priori bound in Proposition 4.5 is sharp in terms of the predicted convergence behavior as we will show now for a test case from [BS18]. Moreover, we will investigate the performance of Algorithm 4.1 also for a test case from [BS18]. On ∂ Ω out \ Γ out and ∂Ω out \Γ out we prescribe homogeneous Neumann boundary conditions and we suppose that Ω out and Ω out do not border the Dirichlet boundary of Ω. For the FE discretization we use a regular mesh with hexahedral elements and a mesh size h = 0.1 in each space direction and a corresponding conforming FE space with linear FE resulting in dim(V h | Ωout ) = 30258, dim(R) = N R = 2172, dim(S) = N S = 2880 for Ω out and V h | Ωout = 55473, N R = 3987, and N S = 5280 for Ω out . 3 We equip the source space S with the L 2 -inner product and the range space R with the energy inner product. Finally, for all results in this subsubsection we computed the statistics over 1000 samples.
Analyzing the convergence behavior of E( T − P R rand k+p T ) on Ω out for a growing number of randomly generated basis functions k and a (fixed) oversampling parameter p = 2 in Fig. 4a we see that until k ≈ 75 the a priori bound reproduces the convergence behavior of E( T − P R rand k+p T ) perfectly. We may thus conclude that    (47), and the a priori error bound of (47) scaled with a constant such that its value for k = 2 equals the one of E( T − P R rand k+p T ) (sc. a priori) for increasing k for and p = 2 for the oversampling domain Ω out (a). Median of the projection error T − P R rand n T for a decreasing target accuracy tol for a varying number of test vectors n t and the minimal and maximal values for n t = 10 on Ω out (b). the a priori bound in (47) seems to be sharp regarding the convergence behavior of E( T − P R rand k+p T ) in the basis size k. We also observe that the a priori bound is rather pessimistic as it overestimates E( T − P R rand k+p T ) by a factor of more than 100; this is mainly due to the square root of the conditions of the inner product matrices.
Regarding the performance of Algorithm 4.1 on Ω out we first observe in Fig. 4b that the actual error T − P R rand n T lies below the target tolerance tol for all 1000 samples for n t = 10; which holds also true for all other considered values of n t . Here, we prescribe ε algofail = 10 −10 and use 3993 as an upper bound for N T . We see in Fig. 4b that increasing the number of test vectors n t from 5 to 10 or from 10 to 20 increases the ratio between the median of the actual error T − P R rand n T and the target accuracy tol significantly -for the former by more than one magnitudewhile an increase from n t = 40 to n t = 80 has hardly any influence; similar results have been obtained in [BS18] for heat conduction and a Helmholtz problem. This can be explained by the scaling of the effectivity of the employed a posteriori error estimator, which is of the order of 1000 for n t = 5 and of the order of 10 for n t ≥ 20. Regarding the choice of n t it seems that for the present test case a value of about 20 is in the sweet spot. We thus infer that for the present test case only very few local solutions in addition to the optimal amount are required, demonstrating that Algorithm 4.1 performs nearly optimally in terms of computational complexity for the current problem. 4.2.3. The general setting P = {μ}. The processes in subsubsection 4.2.1 and 4.2.2 yield for every µ ∈ P the local approximation space R + n (µ) for this specific parameter µ ∈ P. R + n (µ) can also be generated by some other process, where we require that there holds possibly only at high probability and that R + n (µ) is defined as the direct sum of R n (µ), the kernel of the bilinear form, and representations of non-homogeneous Dirichlet boundary conditions and the right-hand side. We abuse notation in this subsubsection by omitting henceforth the remark that the estimate may only hold in a probabilistic sense. The constant C 1 (T H , µ) has to be chosen in such a manner that if one uses the parameter-dependent spaces R + n (µ) to define u n (µ), we have The spectral greedy algorithm as introduced in [SP16] 4 constructs one (quasioptimal) parameter-independent approximation space R N which approximates those parameter-dependent spaces R + n (µ) with a given accuracy on a finite dimensional training set Ξ ⊂ P. In the spectral greedy algorithm we exploit the fact that we expect that the local spaces R + n (µ), and in particular the spectral modes that correspond to the largest eigenvalues, are not affected too much by a variation in the parameter thanks to the expected very rapid decay of the higher eigenfunctions in the interior of Ω out .
The spectral greedy as described in Algorithm 4.2 then proceeds as follows. After the initialization we compute for all µ ∈ Ξ the parameter-dependent spaces R + n (µ) such that we have (48). Note that for a decomposition T H with mutually disjoint interfaces (also called ports), where each Ω m , m = 1, . . . , M has exactly two interfaces and all interfaces have the same geometry we have the following a priori error bound [SP16] for the error between u(µ) and the continuous port-reduced static condensation approximation u n+ (µ) corresponding to the parameter-dependent optimal interface space R + n (µ): Here, the constant C γ,1 (Ω γ , µ) depends only on the subdomains that share γ and not on Ω or on u(µ). Moreover, c 1 (µ) and c 2 (µ) are chosen such that we have c 1 (µ)|||·|||μ ≤ |||·||| µ ≤ c 2 (µ)|||·|||μ for all µ ∈ P and a fixed reference parameterμ ∈ P. Choosing C 1 (T H , µ) = #γc 1 (µ)c 2 (µ) max γ∈T γ H C γ,1 (Ω γ , µ) and λ γ,n+1 (µ) ≤ ε/2 yields a reduced space R + n (µ) that satisfies the requirements stated in the beginning for every µ ∈ Ξ. Although precise estimates for C γ,1 (Ω γ , µ) can be obtained, setting C γ,1 (Ω γ , µ) = 1 yields in general good results as another value would just result in rescaling ε; for further details see [SP16]. After having collected all functions on Γ m,m or Ω m that are essential to obtain a good approximation for all local solutions u loc (µ) of the PDE evaluated on Γ m,m or Ω m , µ ∈ Ξ, we must select a suitable basis from those functions. This is realized in an iterative manner in Lines 5-14. 4 For a generalization to a setting where the discrete parameter set describes different geometries such as a beam with or without a crack we refer to [Sme19].

Algorithm 4.2: spectral greedy [SP16]
Input : train sample Ξ ⊂ P, tolerance ε Output : set of chosen parameters Ξ N , local approximation space R N 1 Initialize N ← dim(ker(a m (·, v))), In each iteration we first identify in Line 7 the reduced space R + n (µ * ) that maximizes the deviation where possible choices of S(R + n (µ)) ⊂ R + n (µ) will be discussed below. Subsequently, we determine in Line 9 the function κ ∈ S(R + n (µ * )) that is worst approximated by the space R N and enhance R N with the span of κ. The spectral greedy algorithm terminates if for all µ ∈ Ξ we have for a constant C 2 (T H , µ), which can in general be chosen equal to one. We emphasize, that both C 1 (T H , µ) and C 2 (T H , µ) do in general only depend on the number of faces or subspaces on which the respective reduced space R N is used and not on the precise decomposition of Ω; see (50). A slight modification of the stopping criterion (51) and a different scaling of ε in the threshold for the a priori error bound in Line 3 allows to prove that after termination of the spectral greedy for a decomposition T H with mutually disjoint interfaces, where each Ω m , m = 1, . . . , M has exactly two interfaces and all interfaces have the same geometry we have [SP16] (52) |||u(µ) − u N (µ)||| µ /|||u(µ)||| µ ≤ ε.
Here, u N (µ) is the continuous port-reduced static condensation approximation corresponding to R N ; R N being the continuous outcome of the spectral greedy.
Choice of the subset S(R + n (µ)) First, we emphasize that in contrast to the standard greedy as introduced in [ properties is implicitly saved in their norms as χ sp j (µ) 2 R = λ j (µ), j = 1, . . . , n. To obtain local approximation spaces R N that yield a (very) good approximation u N (µ) already for moderate N it is therefore desirable that the spectral greedy algorithm selects the lower eigenmodes sooner rather than later during the while-loop. As suggested in [SP16] we thus propose to consider , n + := dim(R + n (µ)) and here and henceforth {χ i (µ)} n+ i=1 denotes the orthonormal basis of R + n (µ). Note that we are therefore considering a weighted norm in R + n (µ). The deviation E(S(R + n (µ)), R N ) can then be computed by solving the eigenvalue problem: Find ( j (µ), σ j (µ)) ∈ (R n+ , R + ) such that and χ k denotes the orthonormal basis of R N . We thus obtain E(S(R + n (µ)), R N ) = σ 1 (µ), for all µ ∈ Ξ, and κ = n+ i=1 1 (µ * )χ i (µ * ) at each iteration. Note that were we to consider the norm · R in (53) the sorting of the spectral basis χ i (µ) of R + n (µ) in terms of approximation properties is neglected in the while loop of Algorithm 4.2; for further explanations see [SP16].
In the online stage we consider E(µ) ≡ 1 in both components, G = (0, 0, 0) T , and prescribe g D,1 = (0, 0, 0) T at Γ 1 and g D,2 = (1, 1, 1) T at Γ 2 . We observe that the Legendre modes perform by far the worst, demonstrating that including information on the solution manifold in the basis construction procedure can significantly improve the approximation behavior. We remark that the Legendre modes will perform even worse in the case of less regular behaviour on the interface, which further justifies the need for problem-specific local approximation spaces in the sense of model reduction. The empirical modes and spectral modes exhibit a comparable convergence until N = 17, but for N > 17 the relative error in the spectral approximation is one order of magnitude smaller than that of the empirical port mode approximation. This can be explained by the fact that thanks to its conception the pairwise training algorithm is able to identify and include the most significant modes, but (in contrast to the spectral greedy algorithm) might have difficulties to detect subtle modes that affect the shape of the function at the interface Γ m,m only slightly. Note that the temporary stagnation of the relative error for N = 7, ..., 17 for the spectral modes is due to the fact that the spectral greedy prepares the interface space for all possible boundary conditions and parameter configurations. Thus, for the boundary conditions considered here some spectral modes, as say a mode related to a twisting (torsion) of the beam, are not needed for the approximation.
For localized model order reduction, however, we are merely interested in localized a posteriori error estimation. To this end, we first present abstract localized lower and upper bounds for the dual norm of a linear functional (see [BEOR17]).
Theorem 5.1 (Localized lower and upper bounds for functionals). Let O i , 1 ≤ i ≤M be a collection of linear subspaces of V h , and let P Oi : Defining the stability constant of this partition modulo V N as Here, f O i denotes the norm of the restriction of f to O i .
When grouping the spaces O i so that in each group, all spaces are orthogonal to each other, J is the number of groups needed. Note that subtracting the projection onto V N in (56) allows subtracting say the mean value of a function or the orthogonal projection onto the rigid body motions, if the respective functions are included in V N . We may thus employ say Poincaré's inequality or Korn's inequality in subdomains that do not lie at Γ D .
Applying both estimates to the residual R(u N (µ); µ) ∈ V h , we obtain from (54) and Theorem 5.1 a robust and efficient, localized error estimate: Corollary 5.2 (Localized residual based a posteriori error estimate). Let the assumptions on the subspace collection O i and the mappings P i from Theorem 5.1 be satisfied. Then, the error estimator ∆ loc (u N (µ)) defined as is robust and efficient, i.e.
Online-offline decomposition of this error estimator can be done by applying the usual strategy for online-offline decomposition used with the standard RB error estimator (see e.g. [HRS16, Sec. 4.2.5] or a numerically more stable approach [BEOR14,CEL14a,SFDE15]) to every dual norm in ∆ loc (u N (µ)).
The a posteriori error estimator for the ArbiLoMod derived in [BEOR17] and the a posteriori error estimator for the scRBE method as suggested in [Sme15] both fit into the framework above as will be detailed below in Examples 5.3 and 5.4. In contrast, for instance the error estimators proposed in [HKP13b,HKP13a] for the scRBE method exploit matrix perturbation analysis at the system level to bound the Euclidean norm of the error between the coefficients of the static condensation solution and the coefficients of the static condensation solution using an RB approximation in the interior. To estimate the error caused by interface reduction in [EP13] a computationally tractable non-conforming approximation to the exact error is employed. To take into account the error due to the intra-element RB approximations ideas from [HKP13b] are used. It can also be noted that the error estimators in [HKP13b,HKP13a,EP13] are only valid under certain assumptions on the accuracy of the RB approximation. In [MRH15] a localized a posteriori error estimator for interface reduction and intra-element RB approximation is presented for the coupled Stokes-Darcy system. The a posteriori error estimator for the CMS method derived in [JBL11] employs the dual norms of residuals and eigenvalues of the eigenproblems used for the construction of the (local) basis functions. The error estimator in [JBL11] is however only partially local as it involves the residual for the port or interface space on the whole interface Γ. For localized a posteriori error estimation in the context of adaptive GMsFEM we refer to [CEL14b,CEL18a,CEL15,CEL17b].
The constant c pu will depend on size of the overlap of the subdomainsΩ i with their neighbors in relation to their diameters.
Moreover, we assume that there is a linear interpolation operator I onto V h such that I is the identity on We then can define mappings which satisfy the assumptions of Theorem 5.1. In case V h comes from a finite element discretization, a possible choice for I is Lagrange interpolation.
If we now ensure that the partition of unity p i is included in V N , we can chooseφ in the definition of c N asφ := p i · |Ω i | −1 Ω i ϕ, which allows us to prove [BEOR17, Proposition 5.7] that c N can be bounded by In this estimate c ovlp := max x∈Ω #{i ∈ Υ E | x ∈ Ω i } is the maximum number of estimator domainsΩ i overlapping in any point x of Ω, and c pc is a Poincaréinequality constant associated withT H . In particular, this result shows that the efficiency of (58) is independent from the number of subdomains in T H , provided that the partition of unity p i is included in V N .
Example 5.4 (scRBE method and interface reduction from [Sme15]). We exemplify the a posteriori error estimator from Corollary 5.2 for the scRBE method, which is equally applicable when considering solely static condensation and no intraelement RB approximations. 7 To simplify notations we define interface spaces ; for the definition of ψ γ k we refer to subsection 3.1. Recall that we then have the following space decomposition of the We may thus uniquely rewrite every ϕ ∈ V h as where ϕ m ∈ V m h;0 and ϕ γ ∈ V γ h , extending ϕ m and ϕ γ by zero. This allows us to define mappings We thus obtain equals the trace of ϕ minus the orthogonal projection on the kernel of the bilinear form; for further details see [BS18,SP16]. Then, we can use [SP16,Lemma B.4] to conclude boundedness of c N and thus (59), the latter corresponding to [Sme15, Proposition 4.2 and Corollary 4.6].
Finally, we shortly discuss how to compute the dual norms of the residuals in (57). The dual norms of the residuals of the intra-element RB approximations can be computed by employing Riesz representations. The dual norms of the residuals in the interface space can be computed by means of conservative fluxes [HEML00], which have been extended to interface reduction in [Sme15]. In detail, we compute the conservative flux H m N (µ) such that where (·, ·) γ denotes a suitable inner product on the interface γ. Note that thanks to our mutual disjoint interface assumption problem (63) decouples and we may compute the conservative flux separately for each interface γ. Moreover, by orthonormalizing the interface basis functions χ γ k defined in subsection 3.1 w.r.t the (·, ·) γ inner product, the computation of H m N (µ) reduces to the assembling of the residual in (63). The computational costs thus scale linearly in (N γ h − N γ ) and N γ . For further details we refer to [Sme15].

5.2.
Local flux reconstruction based error estimation. Following [OS15], we discuss local flux reconstruction based a posteriori error estimation of the full approximation error u(µ) − u N (µ) (that is: the discretization as well as the model reduction error) in the context of non conforming approximations of elliptic multiscale problems such as Example 2.2. An extension to convection-diffusion-reaction problems based on [ESV10] is straightforward. This estimate was introduced in the IP localized non conforming setting of the LRBMS (compare Section 3.2.2).
Recalling the broken Sobolev space and broken gradient operator from Section 3.2.2, the key idea of flux reconstruction based error estimation is to observe that not only the approximate solution u N (µ) is non conforming, but also the approximate diffusive flux −κ(µ)∇ h u N (µ), in the sense that it is not contained in H div (Ω) (i.e., the space of functions in L 2 (Ω) d whose divergence exists in a weak sense and lies in L 2 (Ω)).
We may then obtain computable estimates by comparing these quantities with conforming reconstructions, as detailed further below. 8 The respective reconstructed diffusive flux is locally conservative and is related to the conservative flux reconstruction to compute the dual norm of the residuals in the interface space in Example 5.4.
To begin with, we specify the parameter dependent (semi-)energy norm induced by the bilinear form a for a parameterμ ∈ P, |||·|||μ : H 1 (τ h ) → R, v → |||v|||μ := a(v, v;μ) 1 2 , (by using the broken gradient in the definition of a) and note that we can compare these semi norms for two parameters by means of the affine decomposition of a (compare (68)), with the equivalence constants given by Θ a (µ,μ) := min Qa q=1 Θ q a (µ) Θ q a (μ) −1 and Θ a (µ,μ) := max Qa q=1 Θ q a (µ) Θ q a (μ) −1 , respectively. The first abstract result is the following discretization-agnostic lemma, which leaves the choice of the reconstructions, v and s, open. (We give estimates on the full V h -norm at the end of this subsection.) Lemma 5.5 (Abstract energy norm estimate (Lemma 4.1 in [OS15])). For µ ∈ P, let u(µ) ∈ V denote the weak solution of (1) with the data functions κ and q as in example 2.2. It then holds for arbitrary v N ∈ H 1 (τ h ) andμ ∈ P, that To obtain a fully computable localizable estimate we need to specify the conforming reconstruction of the solution (v in the above lemma) and of the diffusive flux (s in the above lemma). We define both reconstructions w.r.t. the global fine grid τ h and note, that their respective computations can be localized w.r.t. the domain decomposition to allow for offline/online decomposable localized estimates.
We reconstruct the non conforming solution u h (µ) ∈ V h by means of its Oswald interpolant I OS [u h (µ)] ∈ V . We define the corresponding Oswald interpolation operator I OS : V h → V h ∩ V by specifying its values on each Lagrange node ν of τ h : for all nodes which are shared by multiple grid elements, which we collect in τ ν h ⊂ τ h . The definition of the conforming reconstruction of the non conforming diffusive flux −κ(µ)∇ h u h (µ) ∈ L 2 (Ω) d is more involved. Given l ≥ 0, we define the lth order Raviart-Thomas-Nédélec space of vector valued functions by RT N l h (τ h ) := s ∈ H div (Ω) s| t ∈ [P l (t)] d + xP l (t) ∀t ∈ τ h and note that the DoFs of any s h ∈ RT N l h (τ h ) are uniquely defined by specifying the moments of order up to l − 1 of s h | t on all elements t ∈ τ h and the moments of order up to l of s h | σ ·n σ on all faces σ ∈ τ γ h (compare [BF91]). With these preliminaries we define the diffusive flux reconstruction operator R l h : P → [V h → RT N l h (τ h )], given some v h ∈ V h and some µ ∈ P by specifying the DoFs of on all σ ∈ τ γ h and for all ∇r ∈ [P l−1 (t)] d with r ∈ P l (t) on all t ∈ τ h . Given a FOM space V h of polynomial order k ≥ 1, we choose a k − 1st order reconstruction. With this definition, the reconstructed diffusive flux of a given a reduced solution u N (µ) fulfills the following local conservation property, given that the constant function 1 is present in the local reduced spaces V m N : for all Ω m ∈ T H .
When inserting this diffusive flux reconstruction for s in lemma 5.5, this local conservation property is key to obtaining the following estimate.
Theorem 5.6 (Locally computable energy norm a posteriori estimate). Let the domain decomposition T H from Definition 2.5 be such, that the Poincaré-inequality holds on each subdomain Ω m ∈ T H with a constant C m where h m := diam(Ω m ) and where Π m 0 ϕ denotes the mean value of ϕ over Ω m . Let further u(µ) ∈ V be the weak solution of (2) and u N (µ) ∈ V N be the IP localized ROM solution, with 1 ∈ V m N for 1 ≤ m ≤ M . It then holds for arbitraryμ,μ ∈ P that with the a posterior error estimator η(µ;μ;μ) given by We obtain an a posterior error estimate w.r.t. the V h -norm or a full energy norm, |||·||| µ + σ∈τ γ h (·, ·) p σ 1 2 , by noting that u(µ), u(µ) p σ = 0 for a weak solution u(µ) of sufficient regularity.

Basis enrichment and online adaptivity
Model order reduction is usually employed either (i) in the context of realtime decision making and embedded devices, or (ii) in the context of outer-loop applications, such as optimal control, inverse problems or Monte Carlo methods. In (i), one is usually interested in reduced spaces V N of very low dimension to obtain ROMs as small as possible, at the possible expense of very involved offline computations. Here, localized model order reduction may help to reduce the latter, but we can usually not expect the resulting reduced space to be smaller than the one generated using traditional global model order reduction methods. In (ii), however, one is interested in a black-box-like approximation scheme which is queried for a huge amount of parameters, with a somehow "optimal" computational cost (including offline as well as online cost). Here, one may keep high-dimensional data throughout the computational process (offline as well as online), and it is in this context that localized model order reduction techniques may truly outperform other approaches. In the context of PDE constrained optimization this has been investigated e.g. in [OS17,OSS18,WFR19].
The localized a posteriori error estimation as discussed in Section 5 enables adaptive enrichment of the local reduced approximation spaces, whenever the quality of the reduced scheme is estimated to be insufficient -be it due to insufficient training due to lacking computational resources or due to limited knowledge about the range of possible parameters or due to other reasons altogether.
Let us thus assume that an initial (possibly empty) localized reduced approximation space V N is given, compare Section 4. The goal of an adaptive enrichment is to enlarge the local solution spaces with additional modes that reflect non-local influences of the true solution such as channeling effect or singularities. Local adaptive basis enrichment can be employed both offline for the whole parameter range and/or online for a specific chosen parameter. Empirical training followed by offline enrichment is e.g. used in a Greedy manner for the basis construction in ArbiLoMod (cf. Example 5.3) in [BEOR17]. Adaptive enrichment for the GMsFEM is presented in [CEL14b,CEL18a] and online adaptive enrichment in [CEL15,CEL17b]. For the exposition in this section, we restrict to online enrichment as introduced [OS15], i.e. for local enrichment of the basis when a certain parameter is already chosen.
From a birds-eye perspective, we can think of an online adaptive reduced scheme as a p-adaptive FE scheme with problem adapted basis functions, where the local reduced bases are adapted during online enrichment. 9 Thus, we can think of online enrichment in the usual Solve → Estimate → Mark → Refine (SEMR) manner, well known in grid-adaptive discretization schemes. In the Estimate step we employ an a posteriori error estimate η that is localizable w.r.t. the domain decomposition, i.e. η 2 ≤ M m=1 η 2 m , with appropriate local indicators η m . Examples are given in Section 5. As such, most marking strategies from grid-adaptive schemes are applicable, and we give examples in Section 8.1. In this context, refinement is locally done by enriching the local reduced spaces, that is: by adding additional basis functions to the local reduced bases on selected subdomains. We thus presume we are given a parameter µ ∈ P and a reduced solution u N (µ) ∈ V N , the estimated error of which is above a given tolerance.
As an example, we detail the online enrichment procedure used in the context of the LRBMS (compare Section 3.2.2), using the a posteriori error estimation techniques from Section 5.2 . Inspired by domain decomposition as well as numerical multiscale methods, we may then obtain a candidate for the next element of a local reduced basis by solving local corrector problems on a collection T H ⊆ T H of marked subdomains with u N (µ) as boundary values. For each marked subdomain Ω m ∈ T H , we denote byΩ m := Ω m ∈ T H Ω m ∩ Ω m = ∅ an overlapping subdomain and by VΩ m h := v|Ω m v ∈ V h , v| ∂Ωm = 0 the associated restricted FOM space, encoding zero Dirichlet boundary values. We are then looking for a local correction ϕΩ m ∈ VΩ m h , such that where we understand all quantities to be implicitly extended to Ω by zero, if required, and note that ϕΩ m can be computed involving only quantities associated withΩ m . Using this local correction on the overlapping subdomain, we obtain the next 9 We would also like to mention the h-adaptive model order reduction approach from [Car15] which is based on a k-means clustering of the DoFs, but we restrict the exposition here to localization w.r.t. a domain decomposition. Algorithm 6.1: Adaptive online enrichment in the context of the LRBMS.
Input : a marking strategy MARK, an orthonormalization procedure ONB, a localizable offline/online decomposable a posteriori error estimate update all reduced quantities (system matrices, error estimates) w.r.t. the newly added basis elements 10 solve (12) for the reduced solution u N (µ) using the updated quantities 11 return u N (µ) element of the local reduced basis associated with Ω m by an orthonormalization of ϕΩ m + u N (µ) | Ωm with respect to the existing basis on V m N . Given a marking strategy and an orthonormalization procedure, we summarize the adaptive online enrichment used in the context of the LRBMS in Algorithm 6.1.

Computational aspects
In this section we discuss the computational efficiency of localized model order reduction schemes in comparison to standard, non-localized techniques. Imposing a localization constraint on the reduced space naturally yields sub-optimal spaces in the sense of Kolmogorov N -width. However, this is mitigated by the sparse structure of the resulting reduced system matrices. In particular, for problems with large-dimensional parameter domains with localized influence of each parameter component on the solution, we can expect localized ROMs to show comparable or even better online efficiency in comparison to a standard ROM. In addition, localized model order reduction provides more flexibility to balance computational and storage requirements between the offline and online phase and has thus the potential to be optimized with respect to the specific needs. This is particularly favorable for large-scale or multiscale problems, where global snapshot computations are extremely costly or even prohibitive.
In the offline (and enrichment) phase of the localized schemes, only relatively smalldimensional local problems are solved instead of the computation of global solution snapshots. In comparison to a global reduction approach with a parallel solver for snapshot generation (e.g. a domain decomposition scheme), the preparation of the local reduced spaces via training (Section 4) can be performed almost communicationfree, allowing the application of these schemes on parallel compute architecture without fast interconnect such as cloud environments. Via adaptive enrichment of the approximation spaces -based on the solution of local correction problems (Section 6) -smaller and more efficient ROMs can be obtained. In comparison to domain-decomposition methods, where similar correction problems are solved, these correction problems are only solved in regions of the domain where the approximation space is insufficient. Thus, for problems with a localized effect of the parameterization, a significant reduction of the computational effort can be expected in the reduced basis generation process.
In the context of component-based localized model order reduction (e.g. CMS, scRBE, RBHM, RDF) large computational savings can be achieved by the preparation of local approximation spaces (components) w.r.t. arbitrary neighboring components (connected through so-called ports). In addition to parametric changes of the governing equations or computational domain, this allows the (non-parametric) recombination of components in arbitrary new configurations without requiring additional offline computations. 7.1. Online efficiency. In view of Definition 2.8, we can interpret the localized model order reduction methods introduced in Section 3 as standard projection based model reduction methods -such as the reduced basis method -subject to the constraint that the reduced space V N admits a localizing decomposition of the form (11). As such, the usual offline/online decomposition methodology can be applied. To this end, let us assume that the bilinear form a(·, ·; µ) and the source functional f (·; µ) admit affine decompositions for all v, w ∈ V , µ ∈ P with non-parametric bilinear forms a q : V × V → R, functionals f q ∈ V and some parameter functionals Θ a q , Θ f q : P → R. If the given problem is not of the form (68), we can employ empirical interpolation [BMNP04] to compute an approximate affine decomposition.
We begin by computing the reduced approximation space V N using the methods outlined in Section 4. After that, a reduced model is assembled by computing matrix representations A q ∈ R N ×N of a q and vector representations F q ∈ R N of f q w.r.t. a given basis ϕ 1 , . . . , ϕ N of V N , i.e.

(69)
A q ij := a q (v j , w i ), F q i := f q (w i ). After this computationally demanding offline phase, the coordinate representation U N (µ) ∈ R N of the reduced solution u N (µ) of (12) is quickly obtained for arbitrary new parameters µ by solving in the following online phase. The computational effort to determine u N (µ) is of order for the assembly and solution of the dense equation system (70). In particular, we have obtained full offline/online splitting, i.e. the effort to obtain U N (µ) is independent of dim V h . From U N (µ) we can then either reconstruct u N (µ) by linear combination with the reduced basis or evaluate arbitrary linear functionals of u N (µ) by additionally computing vector representations of these functionals in the offline phase. RB methods aim at constructing reduced spaces V N which are near-optimal approximation spaces for the discrete solution manifold {u h (µ) | µ ∈ P} in the sense of Kolmogorov, i.e. it should hold that As d loc N > d N , localized RB methods generally result in larger V N to satisfy a given approximation error tolerance ε. Since we can represent any basis vector of a global RB approximation of (10) w.r.t. the localizing space decomposition (9) as a sum of M tot := M + #T γ H + #T e H + #T v H local vectors, we have the a priori bound d loc M tot ·N < d N . In other words, if we denote by N (N glob ) the number of reduced basis vectors required for a localized (global) RB approximation for given ε and denoting by N loc the maximum dimension of the local RB spaces V m where the constant C T H only depends on the topology of the domain decomposition T H . Whether or not estimate (74) is sharp largely depends on the dependence of the solution u(µ) on the parameter µ. When a change in µ equally affects the solution in all subdomains Ω m , we expect that optimal local RB spaces will be of similar dimension N loc and that N loc ≈ N glob . On the other hand, it may be the case that the influence of µ on u(µ) is weak in many Ω m , in which case N N loc M tot , or that each of the p components of µ ∈ R p affects u(µ) on different subdomains, in which case N loc M tot N glob M tot . Thus, the actual loss in online efficiency due to localization will strongly depend on the type of problem to be solved.
More importantly though, note that the localization of V N results in a change of the structure of the reduced system matrices A q . While these matrices are dense for global RB approximations, localized RB schemes yield A q with a sparse block structure of M tot × M tot blocks of maximum dimension N loc × N loc and a maximum of C cup blocks per row. C cup depends on the specific localization method and on the topology of T H . For instance, for non-conforming methods, C cup − 1 is given by the maximum number of interfaces of a given subdomain Ω m , whereas for the ArbiLoMod with a quadrilateral mesh C cup = 25.
Thus, estimate (74) has to be interpreted in relation to the fact that the computational complexity for solving (70) can be vastly reduced in comparison to (71) by exploiting the structure of the A a . In particular, the costs for assembling (70) can be reduced to O(C T H C cup (N loc ) 2 M ). For the solution of (70) direct or block-preconditioned iterative solvers can be used. For the latter, the computational effort can be expected to increase sub-quadratically in the number of subdomains M . In the scRBE method, the volume degrees of freedom associated with the spaces V m N are eliminated from (70) using static condensation to improve computational efficiency.
7.2. Offline costs and parallelization. While the local RB spaces V m N , V γ N , V e N , V v N can be initialized by decomposing global solution snapshots u h (µ) w.r.t. (9) (see [AHKO12]), the core element of localized RB methods is the construction of local RB spaces from local problems associated with the subdomains Ω m as described in Section 4. This has various computational benefits: First, we can expect a reduction of computational complexity as for most linear solvers we expect a super-linear increase in computational complexity for an increasing dimension of V h , whereas the ratio of the dimensions of V h and the local subspaces in (9) remains constant of order 1/M (for the volume spaces V m h and smaller for the spaces V γ h , V e h , V v h ) as h → 0. Thus, solving O(M N loc ) training problems of size dim V h /M is expected to be faster than solving N global problems. At the same time, we expect N loc to decrease for H → h. Thus smaller subdomains Ω m will generally lead to shorter offline times at the expense of less-optimal spaces V N . In particular, for the non-conforming schemes in Subsection 3.2.2 it is readily seen that for H = h we have V N = V h and that (12) and (10) are equivalent.
Even more important than a potential reduction of complexity is the possibility to choose H small enough such that each local training problem can be solved communication-free on a single compute node without the need for a highperformance interconnect. Also the problem setup and the computation of the reduced system (70) can be performed mostly communication-free: Instead of instantiating a global fine-scale compute mesh, each compute node can generate a local mesh from a geometry definition, solve training problems for a given local RB space and all coupling spaces to obtain corresponding block-entries in A q , F q . Only the local geometry and the resulting reduced-order quantities are communicated (see [BEOR17,Section 8]). This makes localized RB methods attractive for cloudbased environments, where large computational resources can be dynamically made available, but communication speed is limited.
Depending on the problem structure, the use of online enrichment (Section 6) can yield smaller, problem adapted reduced spaces V N . Similar to the training of V N , online enrichment is based on the solution of small independent local problems, that can easily be parallelized. As typically only some fraction of the subspaces of V N undergo enrichment, less computational resources need to be allocated during an online-enrichment phase. It has to be noted, however, that online enrichment leads to a propagation of snapshot data through the computational domain as the value of the current solution u N (µ) at the boundary of the enrichment-problem domain enters the problem definition. Thus, to perform online-enrichment, (boundary values of) reduced basis vectors have to be communicated between compute nodes and the entire reduced basis has to be kept available.

Applications and numerical experiments
8.1. Multiscale problems. We demonstrate the IP localized RB methods from Section 3.2.2 in the context of parametric multiscale problems, such as Example 2.2, with a focus on online adaptivity as in Section 6 (using the a posteriori error estimate from Section 5.2), rather than offline training. These experiments were first published in the context of the online adaptive LRBMS in [OS15]. We consider a multiplicative splitting of the parameter dependency and the multiscale nature of the data functions, in the sense that κ(µ) := λ(µ)κ ε , with a parametric total mobility λ : P → L ∞ (Ω) and a highly heterogeneous permeability field κ ε ∈ L ∞ (Ω) d×d . To be more precise, we consider (2)  On each t ∈ τ h , κ ε | t is the corresponding 0th entry of the permeability tensor used in the first model of the 10th SPE Comparative Solution Project (which is given by 100 × 20 constant tensors, see [SPE01]) and λ c models a channel, as depicted in Fig.  6, top left.
The right hand side f models a strong source in the middle left of the domain and two sinks in the top and right middle of the domain, as is visible in the structure of the solutions (see Fig. 6, third row). The role of the parameter µ is to toggle the existence of the channel λ c . Thus λ(1)κ ε = κ ε while µ = 0.1 models the removal of a large conductivity region near the center of the domain (see the second row in Fig. 6). This missing channel has a visible impact on the structure of the pressure distribution as well as the reconstructed velocities, as we observe in the last two rows of Fig. 6. With a contrast of 10 6 in the diffusion tensor and an ε of about |Ω|/2,000 this setup is a challenging heterogeneous multi-scale problem.
We used several software packages for this numerical experiment and refer to [OS15] for a full list and instructions on how to reproduce these results. We would like to mention that all grid-related structures (such as data functions, operators, functionals, products, norms) were implemented in a DUNE-based C++ discretization (which is by now contained in the DUNE extension modules 10 and the generic discretization toolbox dune-gdt 11 ), while we used pyMOR [MRS16] for everything related to model reduction (such as Gram-Schmidt, Greedy). We consider a domain decomposition of |T H | = 25 × 5 squares, each refined such that the full global grid would consist of |τ h | = 1,014,000 elements. For the IP localized FOM, following Section 3.2.2, we choose on each subdomain Ω m ∈ T H the DG space (1st order), product and bilinear form from example 3.2. For error estimation, we employed the flux reconstruction ansatz from Section 5.2 using a zero order diffusive flux reconstruction (compare Theorem 5.6).
The sole purpose of these experiments is to demonstrate the capabilities of localized RB methods regarding online enrichment. We thus initialize the local reduced spaces V m N on each subdomain a priori by orthonormalized Lagrangian shape functions of order up to one, thus obtaining a reduced space with poor approximation properties (comparable to a standard DG space w.r.t. the domain decomposition). Since we employ the a posteriori error estimate η on the full approximation error (including the discretization as well as the model reduction error) from Theorem 5.6, and since we omit grid-refinement in these experiments, the estimated discretization error over all parameters of 1.66 is a lower bound for the overall approximation error, and we thus choose a tolerance of ∆ online = 2 for the online enrichment in Algorithm 6.1.
We compare two different strategies, corresponding to the two plots in Fig. 7. In both cases, we simulate an outer-loop application in the online part by randomly choosing ten parameters P online ⊂ P which are subsequently processed. For each parameter, the local reduced spaces are enriched according to Algorithm 6.1 and the respective marking strategy, until the estimated error is below the specified tolerance. Note that the evaluation of the localizable a posteriori error estimate can be fully offline/online decomposed and that after each enrichment only information from a subdomain and its neighbors are required to locally update the offline/online decomposed data.
In the first experiment, we use a uniform marking strategy, which results in an unconditional enrichment on each subdomain (comparable to domain decomposition methods). As we observe in Fig. 7 (top), however, it takes 129 enrichment steps to lower the estimated error below the desired tolerance for the first online parameter µ 0 . After this extensive enrichment it takes 12 steps for µ 1 and none or one enrichment steps to reach the desired tolerance for the other online parameters.
The resulting coarse reduced space is of size 10,749 (with an average of 86 basis functions per subdomain), which is clearly not optimal. Although each subdomain was marked for enrichment, the sizes of the final local reduced bases differ since the local Gram Schmidt basis extension may reject updates (if the added basis function is locally not linearly independent). As we observe in Fig. 8 (left) this is indeed the case with local basis sizes ranging between 24 and 148. Obviously, a straightforward domain decomposition ansatz without suitable training is not feasible for this setup. This is not surprising since the data functions exhibit strong multiscale features and non-local high-conductivity channels connecting domain boundaries, see Fig. 6.
To remedy the situation we allow for two global snapshots during the offline phase (for parameters µ ∈ {0.1, 1}) and use an adaptive marking strategy which combines uniform marking, Dörfler marking and age-based marking (see the caption of Fig.  7) in the online phase. This strategy employs uniform marking until a saturation condition is reached, and afterwards uses a Dörfler marking combined with a marking based on counting how often a subdomain has not been marked. With two global solution snapshots incorporated in the basis the situation improves significantly, as we observe in Fig. 7 (bottom left). In total we observe only two enrichment steps with uniform marking (see the first two steps for µ 0 ), which indicates that further offline training would be beneficial. The number of elements marked range between 11 and 110 (over all online parameters and all but the first two enrichment steps) with a mean of 29 and a median of 22. Of these marked elements only once have 87 out of 110 elements been marked due to their age (see the last step for µ 1 ). Overall we could reach a significantly lower overall basis size than in the previous setup (1,375 vs. 10,749) and the sizes of the final local bases range between only nine and 20 (compared to 24 to 148 above). We also observe in Fig. 8 (right) that the spatial distribution of the basis sizes follows the spatial structure of the data functions (compare Fig. 6), which nicely shows the localization qualities of our error estimator. 8.2. Fluid dynamics. Flow simulations in pipelined channels have a growing interest in many biological and industrial applications. The localized model order reduction approaches presented in this chapter are suitable for the study of internal flows in hierarchical parametrized geometries. In particular, the non-conforming approach introduced in Section 3.2 has applications in the analysis of the blood flow in specific compartments of the circulatory system that can be represented as a combination of few deformed vessels from a reference one.
We want to solve the Stokes equation defined in (3), with δ = 0, in a computational domain Ω composed by two stenosed blocks Ω µ1 and Ω µ2 (Fig. 9), by imposing  non-homogeneous BCs σ in n = [0, 5] T in the inlet surface (x 1 = 10), non-homogeneous BCs σ in n = [0, −1] T in the outlet surface (x 1 = 0) and homogeneous Dirichlet BC on the remaining boundaries of the domain. Here, the Taylor-Hood Finite Element Method has been used to compute the basis functions, P 2 elements for velocity and supremizer, P 1 for pressure, respectively and consequently P 1 (Γ m,m ) for the Lagrange multipliers space. Fig. 11 shows the distribution of the parameter values selected by the greedy algorithm, by applying the offline stage of the reduced basis method to the single stenosis block. By taking into account that the range [−5, 5] is not admitted, we can see that the higher concentration of values is in the intervals [−10, −5] and [5, 10] in correspondence to larger deformation of the pipe.
The geometry of a single stenosis is obtained by the deformation of a reference pipe through a parameter that represents the contraction in the middle of the pipe. The deformed domain Ω µ is mapped from the straight reference pipeΩ of length L = 5 and radius r = 1 through the following coordinate transformation T µ :Ω → Ω µ such  as x = T µ (x) and x 1 =x 1 +x 1 µ (cos( 2πx3 L ) − 1), x 2 =x 2 +x 2 µ (cos( 2πx3 L ) − 1), x 3 =x 3 . The range of the parameter µ is [−20, −5] ∪ [5, 20], Fig. 10 shows the reference pipe and some representative deformations of the geometry. In order to compute the basis functions, we consider a parametrized Stokes problem for each subdomain. For the first subdomain, we compute the reduced basis imposing zero Dirichlet condition on the wall, Neumann boundary conditions given by imposing σ n = σ · n = ν ∂u ∂n − pn to be σ in n = [0, 5] T on Γ in and σ out n = 0 on the internal interface. For the second subdomain, we compute the reduced basis imposing zero Dirichlet condition on the wall, Neumann boundary conditions imposing σ in n = 0 on the internal interface and σ out n = [0, −1] T on the outflow interface Γ out . Moreover, we enrich the local RB spaces by a coarse finite element solution of the problem computed in the global domain. This strategy ensures not only the continuity of the velocity, but also the one of the normal stress along the internal interface. For this reason this method is called reduced basis hybrid method. Coarse and fine grids have been chosen in order to deal with respectively 155 and 2714 nodes in a single block domain. Fig. 12 shows a representative flow solution in Ω, found with the reduced basis hybrid method, to be compared with the finite element solution. The same comparison, regarding the pressure solutions, is shown in Fig.  13. 9. Further perspectives 9.1. Parabolic problems. Most of the techniques presented in this chapter so far can be extended or even directly applied to parabolic problems. For instance, local approximation spaces that are optimal in the sense of Kolmogorov are proposed in [Sch19] and the LRBMS for parabolic problems is presented in [ORS17,OR17]. To facilitate an adaptive construction of the local reduced space or online-adaptivity, a suitable, localized a posteriori error estimator is key. Therefore, we present in this subsection an abstract framework for a posteriori error estimation for approximations of scalar parabolic evolution equations, based on elliptic reconstruction techniques. For further reading and the application to localized model reduction we refer to [GLV11,ORS17].
Depending on the error we want to quantify, the space V in (75) can be either a analytical function space as in (1) or an already discretized function space V h . We drop the parameter dependency in this section to simplify the notation.
Definition 9.2 (Approximations of the parabolic problem). Let V ⊆ H be a finite dimensional approximation space for V , not necessarily contained in V . Potential candidates for V are conforming or non-conforming localized model reduction spaces V N as discussed above, but also finite element or finite volume spaces fit into this setting. Denote by (·, ·), · the H-inner product and the norm induced by it.
Let f ∈ H, and let a h : (V + V ) × (V + V ) → R be a discrete bilinear form which coincides with a on V × V and is thus continuous and coercive on V . Let further |||·||| be a norm over V + V , which coincides with the square root of the symmetric part of a h over V .
From the definition it is clear thatṽ is the V -Galerkin approximation of the elliptic reconstruction R ell (ṽ).
It is straightforward to modify the estimate in Theorem 9.4 for semi-discrete solutionsũ(t) to take the time discretization error into account: Corollary 9.5. Letũ ∈ L 2 (0, T end , V ), ∂ tũ ∈ L 2 (0, T end , V ) be an arbitrary discrete approximation of u(t), not necessarily satisfying (76). Let R T [p](t) ∈ V denote the V -Riesz representative w.r.t. the H-inner product of the time-stepping residual of u(t), i.e.
In the context of localized model order reduction empirical interpolation has been employed in, e.g., [CEGG14,PY15,OR17]. Based on the concept of empirical operator interpolation from [DHO12] localization strategies can be employed as follows. To present the main ideas, let us assume the simple situation that , v m h ; µ), with a h u h (µ), ·; µ ∈ (V h ) . The strategy will then rely on an empirical operator interpolation of the local volume operators a m h (u m h (µ), ·; µ) ∈ (V m h ) and will thus only involve localized computations in the construction of the interpolation operator. As an example, the interpolation of the local volume operator will be of the form