# Error control for the localized reduced basis multi-scale method with adaptive on-line enrichment

\SetWatermarkTextPREPRINT \SetWatermarkLightness0.85 \SetWatermarkScale4 \setkomafonttitle \setkomafontsection \setkomafontsubsection

Abstract. In this contribution we consider localized, robust and efficient a-posteriori error estimation of the localized reduced basis multi-scale (LRBMS) method for parametric elliptic problems with possibly heterogeneous diffusion coefficient. The numerical treatment of such parametric multi-scale problems are characterized by a high computational complexity, arising from the multi-scale character of the underlying differential equation and the additional parameter dependence. The LRBMS method can be seen as a combination of numerical multi-scale methods and model reduction using reduced basis (RB) methods to efficiently reduce the computational complexity with respect to the multi-scale as well as the parametric aspect of the problem, simultaneously. In contrast to the classical residual based error estimators currently used in RB methods, we are considering error estimators that are based on conservative flux reconstruction and provide an efficient and rigorous bound on the full error with respect to the weak solution. In addition, the resulting error estimator is localized and can thus be used in the on-line phase to adaptively enrich the solution space locally where needed. The resulting certified LRBMS method with adaptive on-line enrichment thus guarantees the quality of the reduced solution during the on-line phase, given any (possibly insufficient) reduced basis that was generated during the offline phase. Numerical experiments are given to demonstrate the applicability of the resulting algorithm with online enrichment to single phase flow in heterogeneous media.

Key words. model reduction, multi-scale, error analysis, flux reconstruction.

AMS subject classifications. 65G99, 65N55, 65N15, 35J20, 65N30, 76S05.

## 1 Introduction

We are interested in efficient and reliable numerical approximations of elliptic parametric multi-scale problems. Such problems consist of finding , such that

(1) |

in a suitable space for parameters , for , where the data functions involved may depend on an a-priori given multi-scale parameter (described in detail in Section 2). An approximation of is usually obtained by a discretization of (1) resulting from a Galerkin projection onto a discrete space , associated with a triangulation of . Since should resolve all features of (1) associated with the fine scale , solving parametric heterogeneous multi-scale problems accurately can be challenging and computationally costly, in particular for strongly varying scales and parameter ranges (see the references in [44]). Two traditional approaches exist to reduce the computational complexity of the discrete problem: numerical multi-scale methods and model order reduction techniques. Numerical multi-scale methods reduce the complexity of multi-scale problems with respect to for a fixed , while model order reduction techniques reduce the complexity of parametric problems with respect to for moderate scales (see [44] for an overview). In general, numerical multi-scale methods capture the macroscopic behavior of the solution in a coarse approximation space , e.g., associated with a coarse triangulation of , and recover the microscopic behavior of the solution by local fine-scale corrections. Inserting this additive decomposition into (1) yields a coupled system of a fine- and a coarse-scale variational problem. By appropriately selecting trial and test spaces and defining the localization operators to decouple this system, a variety of numerical multi-scale methods can be recovered, e.g., the multi-scale finite element method (MsFEM) [22], the variational multi-scale method [43, 37], the multi-scale finite volume method [30] and the heterogeneous multi-scale method (HMM) [5, 31], just to name a few. Model order reduction using reduced basis (RB) methods, on the other hand, is based on the idea to introduce a reduced space , spanned by discrete solutions for a limited number of parameters . These training parameters are iteratively selected by an adaptive Greedy procedure (see e.g. [56] and the references therein). Depending on the choice of the training parameters and the nature of the problem, is expected to be of a significantly smaller dimension than . Additionally, if allows for an affine decomposition with respect to , its components can be projected onto , which can then be used to effectively split the computation into an off-line and on-line part. In the off-line phase all parameter-independent quantities are precomputed, such that the on-line phase’s complexity only depends on the dimension of . The idea of the recently presented localized reduced basis multi-scale (LRBMS) approach (see [36, 8]) is to combine numerical multi-scale and RB methods and to generate a local reduced space for each coarse element of , given a tensor product type decomposition of the fine approximation space, . The coarse reduced space is then given as , resulting in a multiplicative decomposition of the solution into , where the RB functions capture the microscopic behavior of the solution associated with the fine scale and the coefficient functions only vary on the coarse partition .

Other model reduction approaches for parametric (but not multi-scale) problems that incorporate localization ideas and concepts from domain decomposition (DD) methods are the reduced basis element method [39, 38], the reduced basis hybrid method [33, 34] and the port reduced static condensation reduced basis element method [52]. While the idea of the former is to share one reduced basis on all subdomains the idea of the latter two is to generate one reduced basis for each class of subdomains which are then coupled appropriately. There also exist several approaches to use model reduction techniques for homogenization problems (see [14]) and problems with multiple scales, such as the reduced basis finite element HMM [3, 4]. For the case of no scale-separation there exist the generalized MsFEM method [21], which incorporates model reduction ideas, and most recently a work combining the reduced basis method with localized orthogonal decomposition (see [6]).

It is vital for an efficient and reliable use of RB as well as LRBMS methods to have access to an estimate on the model reduction error. Such an estimate is used to drive the adaptive Greedy basis generation during the off-line phase of the computation and to ensure the quality of the reduced solution during the on-line phase. It is usually given by a residual based estimator involving the stability constant and the residual in a dual norm. It was shown in [8] that such an estimator can be successfully applied in the context of the LRBMS, but it was also pointed out that an estimator relying on global information might not be computationally feasible since too much work is required in the off-line part of the computation.

The novelty of this contribution lies in a completely different approach to error estimation – at least in the context of RB methods. We make use of the ansatz of local error estimation presented in [26] which measures the error by a conforming reconstruction of the physical quantities involved, specifically the diffusive flux . This kind of local error estimation was proven to be very successful in the context of multi-scale problems and very robust with respect to . We show in this work how we can transfer those ideas to the framework of the LRBMS to obtain an estimate of the error . We would like to point out that we are able to estimate the error against the weak solution in a parameter dependent energy norm while traditional RB approaches only allow to estimate the model reduction error in a parameter independent norm and only against the discrete solution. In principal, this approach is able to turn the LRBMS method into a full multi-scale approximation scheme, while traditional RB methods can only be seen as a model reduction technique. We would also like to point out that, to the best of our knowledge, this approach (first published in [45]) is the first one to make use of local error information in the context of RB methods.

## 2 Problem formulation and discretization

We consider linear elliptic problems in a bounded connected domain , for , with polygonal boundary and a multiplicative splitting of the influences of the multi-scale parameter and the parameter . An example is given by the problem of finding a global pressure for a set of admissible parameters , such that

(2) |

in a weak sense in , where denotes the usual Sobolev space of weakly differentiable functions and its elements which vanish on the boundary in the sense of traces. Problems of this kind arise for instance in the context of the global pressure formulation of two-phase flow in porous media. Using an IMPES discretization scheme of the pressure/saturation system (see [17] and the references therein), an equation of the form (2) has to be solved in each time step for a different parameter . In that context denotes the scalar total mobility and denotes a possibly complex heterogeneous permeability tensor; external forces are collected in and the parameter models the influence of the global saturation (see Definition 2.1 for details on , and ). Note that more complex boundary conditions and additional parameter dependencies of the boundary values as well as the right hand side (modeling parameter dependent wells, for instance) are possible, but do not lie within the scope of this work.

Triangulation. We require two nested partitions of , a coarse one, , and a fine one, . Let be a simplicial triangulation of with elements . In the context of multi-scale problems we call a fine triangulation if it resolves all features of the quantities involved in (2), specifically if is constant for all . This assumption is a natural one in the context of two-phase flow problems where the permeability field is usually given by piecewise constant measurement data. For simplicity we require to fulfill the requirements stated in [26, Sect. 2.1], namely shape-regularity and the absence of hanging nodes; an extension to more general triangulations is possible analogously to [26, A.1]. We only require the coarse elements to be shaped such that a local Poincaré inequality in is fulfilled (see the requirements of Theorem 4.2). We collect all fine faces in , all coarse faces in and denote by and the neighbors of and , respectively and by the diameter of any element of the sets , , or . We collect in the fine elements of that cover the coarse element and in all faces that cover the set , e.g. by the faces of a fine element , by the faces that cover a coarse face and so forth; the same notation is used for coarse faces . In addition we denote the set of all boundary faces by and the set of all inner faces, that share two elements, by , such that and . We also denote the set of fine faces which lie on the boundary of any coarse element by and by the set of fine faces which lie in the interior of the coarse element. Finally, we assign a unit normal to each inner face , pointing from to , and also denote the unit outward normal to by for a boundary face , for .

### 2.1 The continuous problem

For our analysis we define the broken Sobolev space by as the most general space for the weak, the discrete and the reduced solution defined below. In the same manner we denote the local broken Sobolev spaces for all coarse elements and denote by the broken gradient operator which is locally defined by for all and . Now given , and as stated in Definition 2.1 we define the parametric bilinear form , , and the linear form by and , respectively, and their local counterparts , , and for all , and by

and |

respectively.

###### Definition 2.1 (Weak solution).

Let be bounded, be strictly positive for all and symmetric and uniformly positive definite, such that is bounded from below (away from 0) and above for all . We define the weak solution for a parameter , such that

(3) |

Note that, since is continuous and coercive for all (due to the assumptions on ) and since is bounded, there exists a unique solution of (3) due to the Lax-Milgram Theorem. Given these requirements we denote by and the smallest and largest eigenvalue of , respectively, for any and and additionally define and .

A note on parameters. In addition to the assumptions we posed on above we also demand it to be affinely decomposable with respect to , i.e there exist strictly positive coefficients and nonparametric components , for , such that . We can then compare for two parameters by , where and denote the positive equivalence constants. This assumption on the data function is a common one in the context of RB methods and covers a wide range of physical problems. If does not exhibit such a decomposition one can replace by an arbitrary close approximation using Empirical Interpolation techniques (see [10, 20]) which does not impact our analysis. All quantities that linearly depend on inherit the above affine decomposition in a straightforward way. In particular there exist component bilinear forms such that and (and their discrete counterparts introduced further below) are also affinely decomposable. Since we would like to estimate the error in a problem dependent norm we introduce a parameter dependent energy (semi) norm , , defined by with the local semi-norm defined by , for all , and ; the local semi-norm on the fine triangulation, for any , is defined analogously to . Note that is a norm only on . We can compare these norms for any two parameters using the above decomposition of (the same holds true for the local semi-norms):

(4) |

### 2.2 The discretization

We discretize (3) by allowing for a suitable discretization of at least first order inside each coarse element and by coupling those with a Symmetric weighted Interior-Penalty discontinuous Galerkin (SWIPDG) discretization along the coarse faces of . This ansatz can be either interpreted as an extension of the SWIPDG discretization on the coarse partition , where we further refine each coarse element and introduce an additional local discretization, or it can be interpreted as a domain-decomposition approach, where we use local discretizations, defined on subdomains given by the coarse partition, which are then coupled by the SWIPDG fluxes. In view of the latter, this ansatz shows some similarities to [16] but allows for a wider range of local discretizations. A similar ansatz for a multi-numerics discretization using a different coupling strategy was independently developed and recently introduced in [46]. We will present two particular choices for the local discretizations and continue with the definition of the overall DG discretization.

Local discretizations. The main idea is to approximate the local bilinear forms , which are defined on the local subdomain triangulations , by local discrete bilinear forms discretizing (2) on with homogeneous Neumann boundary values. We will additionally choose local discrete ansatz spaces , with polynomial order , to complete the definition of the local discretizations. The first natural choice for the local discretization is to use a standard continuous Galerkin (CG) discretization, which we obtain by setting and with

where denotes the set of polynomials on with total degree at most . Another choice is to use a discontinuous space , given by | ||||

and to choose from a family of DG discretizations. Therefore we introduce the technicalities needed to state a common framework for the non symmetric, the incomplete, the symmetric and the symmetric weighted interior-penalty (IP) DG discretization (henceforth denoted by NIPDG, IIPDG, SIPDG and SWIPDG, respectively, see [27] and the references therein), following [26, Sect. 2.3].

For a function that is double-valued on interior faces we denote its jump on an inner face by with , recalling that for . We also assign weights to each inner face, such that , and denote the weighted average of by . On a boundary face we set , , and . With these definitions we define the local bilinear form , for by

(5) |

on , with its coupling and penalty parts and , respectively, defined by

and |

for all , all and all . The parametric positive penalty function is given by , where denotes a user-dependent parameter and the locally adaptive weight is given by for an interior face and by on a boundary face , respectively, with . Using the weights we obtain the NIPDG bilinear form for , the IIPDG bilinear form for and the SIPDG bilinear form for . We obtain the SWIPDG bilinear form for by using locally adaptive weights and . From now on we assume that is of the form (5), since this is the most general case and also covers a CG discretization, where all face terms vanish due to the nature of .

Global coupling. Now given suitable local discretizations on the coarse elements we couple those along the coarse faces using a SWIPDG discretization again and define the bilinear form , , by

(6) | ||||

where we use the SWIPDG variants of to define the coupling bilinear form , , by | ||||

(7) |

for all , all and all . Accordingly, we define the DG space for by

with either being the local CG space or the local DG space .

###### Definition 2.2 (DG solution).

Let be large enough, such that is continuous and coercive for all . We define the DG solution for a parameter , such that

(8) |

As always with DG methods, coercivity is to be understood with respect to a DG norm on (i.e., given by the semi norm combined with a DG jump norm) and there exists a unique solution of (8) due to the Lax-Milgram Theorem, if is chosen accordingly.

###### Remark 2.3 (Properties of the discretization).

Depending on the choice of the coarse partition and the local discretization we can recover several discretizations from (6). Choosing and yields a standard CG discretization (except for the treatment of the boundary values), for instance. Choosing , on the other hand, results in the standard SWIPDG discretization proposed in [27, 26]. Note that the local discretizations as well as the polynomial degree do not have to coincide on each coarse element (or even inside a coarse element when using a local DG space). It is thus possible to balance the computational effort by choosing local CG or -adaptive DG discretizations. This puts our discretization close to the multi-numerics discretization proposed in [46], where the latter allows for an even wider range of local discretizations while coupling along the coarse faces using Mortar methods. Concerning the choice of the user dependent penalty factor, we found an automated choice of depending on the polynomial degree , as proposed in [23], to work very well.

## 3 Model reduction

Disregarding the multi-scale parameter for the moment, model reduction using reduced basis (RB) methods is a well established technique to reduce the computational complexity of (8) with respect to the parameter . It has been successfully applied in multi-query contexts, where (8) has to be solved for many parameters (e.g. in the context of optimization), and in real-time contexts, where a solution of (8) (or a quantity depending on it) has to be available for some as fast as possible (see [8, 20, 44] and the references therein). The general idea of RB methods is to span a reduced space by post-processed solutions of (8) for a selection of parameters and to apply a Galerkin projection of (8) onto this reduced space. Together with the assumption we posed on the data function in section 2 this allows for the well-known off-line/on-line decomposition of the RB method, where all quantities depending on (in particular high-dimensional evaluations of the bilinear form) can be precomputed and stored in so called Gramian-matrices during a possibly computationally expensive off-line phase. In the on-line phase of the simulation only low dimensional quantities depending on the parameter dependency of the data functions and the dimension of need to be computed which can usually be done in real time and even on low end devices such as smartphones or embedded devices.

It has been shown in [8], however, that the computational complexity of the off-line phase can become unbearable if the RB method is applied to parametric multi-scale problems, such as (8) for small scales . There are two main drawbacks of classical RB methods in this context: all high-dimensional quantities (solution snapshots, operator matrices, …) usually depend on the fine triangulation , the size of which scales with . Correspondingly all evaluations involving these quantities (scalar products, operator inversions, …) become increasingly expensive for small scales ; the generation of the reduced space usually requires the discrete problem (8) to be solved for several parameters . In the context of multi-scale problems, however, it is usually only feasible to carry out very few global solutions of the original problem on the fine triangulation, if at all.

The first of these shortcomings has been successfully addressed by the localized reduced basis multi-scale (LRBMS) method introduced in [36, 8], which takes advantage of the coarse partition and the underlying discretization to form local reduced spaces on each coarse element instead of the usual global RB space. This allows to balance the computational effort of the off-line versus the on-line phase by selecting an appropriate size of the coarse partition depending on the multi-scale parameter and the real-time or multi-query context of the application. But it was also shown in [8], that the standard residual based estimator which is usually used in the context of RB methods does not scale well in the context of multi-scale problems. We thus finalize the definition of the LRBMS method in this contribution with an efficiently off-line/on-line decomposable a-posteriori error estimator (see the next section). The second shortcoming of classical RB methods is addressed in section 5, where we propose an on-line enrichment extension for the LRBMS method based on the presented error estimator, which finally turns the LRBMS method into a fully capable multi-scale method.

The reduced problem. Let us assume that we are given a local reduced space on each coarse element . Those are usually, but not necessarily, spanned by solutions of (8) and of low dimension, i.e., . We then define the coarse reduced space by

and obtain the reduced solution by a Galerkin projection of (8) onto .

###### Definition 3.1 (Reduced solution).

We define the reduced solution for a parameter , such that

(9) |

Off-line/on-line decomposition. The well-known off-line/on-line decomposition of RB methods relies on an affine parameter dependence of the data functions (see section 2) that carries over to the discrete bilinear form, i.e., there exist nonparametric component bilinear forms , for , such that , with being the coefficients of the corresponding decomposition of . In the standard RB framework, given a basis of , with , one would precompute dense Gramian matrices , given by during the off-line phase. In the on-line phase, given a parameter , one would obtain the reduced system matrix by a low dimensional summation , which only involves evaluations of the scalar coefficient functions for the current parameter. The resulting dense reduced system is of size and does not depend on the dimension of . It can thus be solved with a complexity of in the on-line phase. Denoting by and (with ) the number of Degrees of Freedom of the global and the local spaces, respectively, the reduction of during the off-line phase, however, would be of a computational complexity of , which can become too expensive in the context of small scales .

The flexible framework of the LRBMS method, however, allows for a more local approach. Since the affine decomposition of the data function also carries over to the local bilinear forms and the coupling bilinear forms , we define local Gramian matrices on each coarse element , given by

where denotes a basis of with . In addition, we define coupling Gramian matrices for all coarse elements and their neighbors by

In the same manner, we define local vectors , given by , for all . The global Gramian matrices and the global vector , with , are then given by a standard DG mapping (with respect to the coarse triangulation) of their local counterparts. The complexity to directly solve the reduced system in the on-line phase, , is by definition larger than for standard RB methods. The reduction of , however, can be carried out with a complexity of roughly , which can be computed significantly faster than in the standard RB setting described above (since local computations can be carried out in parallel), depending on the choice of .

###### Remark 3.2 (Properties of the reduced system).

The reduced system matrix of the LRBMS method is obviously larger than the one we would obtain for standard RB methods (given that it scales with ). On the other hand we obtain a sparse matrix of dense blocks (stemming from the local and coupling blocks and the DG mapping) that is usually much smaller than the system matrix of the high-dimensional problem. Thus, if is large, the reduced system can also be solved using sparse iterative solvers with a complexity of .

## 4 Error analysis

Our error analysis is a generalization of the ansatz presented in [26] to provide an estimator for our DG approximation (8) as well as for our LRBMS approximation (9). The main idea of the error estimator presented in [57, 26] is to observe that the approximate DG diffusive flux is nonconforming while its exact counterpart belongs to , which denotes the space of vector valued functions the divergence of which lies in . The idea of [57, 26] is to reconstruct the discrete diffusive flux in a conforming Raviart-Thomas-Nédélec space and compare it to the nonconforming one. Their error analysis relies on a local conservation property of the reconstructed flux on the fine triangulation to prove estimates local to the fine triangulation.

We transfer this concept to the DG discretization defined in section 2.2 and prove estimates local to the coarse triangulation that are valid for the DG approximation as well as for our LRBMS approximation. We obtain mild requirements for the coarse triangulation and the local approximation spaces, namely that a local Poincaré inequality holds on each coarse element and that the constant function is present in the local approximation spaces. The latter is obvious for traditional discretizations and can be easily achieved for the LRBMS approximation. The estimates are fully off-line/on-line decomposable and can thus be efficiently used for model reduction in the context of the LRBMS. Preliminary results were published in [45].

We begin by stating an abstract energy norm estimate (see [26, Lemma 4.1]) that splits the difference between the weak solution of (3) and any function into two contributions. This abstract estimate does not depend on the discretization and thus leaves the choice of and open. Note that we formulate the following Lemma with different parameters for the energy norm and the weak solution. The price we have to pay for this flexibility is the additional constants involving and , that vanish if and coincide. In the following we denote the product over a space , for , by and omit if ; the same holds for the induced norm .

###### Lemma 4.1 (Abstract energy norm estimate).

Let be the weak solution of (3) for and let and be arbitrary. Then

(10) | ||||

###### Proof.

We mainly follow the proof of [26, Lemma 4.1] while accounting for the parameter dependency of the energy norm and the weak solution. It holds for arbitrary , and , that

(11) |

(see [57, Lemma 7.1]) and for the weak solution of (3), that

(12) |

for all and all , where we used the definition of in the first equality and the fact that due to Green’s Theorem and in the second one. Inserting (12) into (11) with and using the norm equivalence (4) then yields the first inequality in (10).

To obtain the second inequality we choose and in the right hand side of (10) which eliminates the two infimums and leaves us with two terms yet to be estimated arising inside the supremum. Using Green’s Theorem and the definition of we observe the vanishing of the first term. We estimate the second term as

where we used the Cauchy-Schwarz inequality and the definition of the energy norm. We finally obtain the second inequality of (10) from the bound above by observing that the supremum vanishes (due to ) and by using the norm equivalence (4) again. ∎

The following theorem states the main localization result and gives an indication on how to proceed with the choice of : it allows us to localize the estimate of the above Lemma, if fulfills a local conservation property. It is still an abstract estimate in the sense that it does not use any information of the discretization and does not yet fully prescribe and .

###### Theorem 4.2 (Locally computable abstract energy norm estimate).

Let be the weak solution of (3) for , let and be arbitrary, let fulfill the local conservation property and let denote the constant from the Poincaré inequality for all on all , where denotes the -orthogonal projection onto for and . It then holds for arbitrary , that

with the abstract global estimator defined as

and the local nonconformity estimator defined as , the local residual estimator defined as and the local diffusive flux estimator defined as for all coarse elements , where .