## Abstract

A general framework to combine numerical homogenization and reduced-order modelling techniques for partial differential equations (PDEs) with multiple scales is described. Numerical homogenization methods are usually efficient to approximate the effective solution of PDEs with multiple scales. However, classical numerical homogenization techniques require the numerical solution of a large number of so-called microproblems to approximate the effective data at selected grid points of the computational domain. Such computations become particularly expensive for high-dimensional, time-dependent or nonlinear problems. In this paper, we explain how numerical homogenization method can benefit from reduced-order modelling techniques that allow one to identify offline and online computational procedures. The effective data are only computed accurately at a carefully selected number of grid points (offline stage) appropriately ‘interpolated’ in the online stage resulting in an online cost comparable to that of a single-scale solver. The methodology is presented for a class of PDEs with multiple scales, including elliptic, parabolic, wave and nonlinear problems. Numerical examples, including wave propagation in inhomogeneous media and solute transport in unsaturated porous media, illustrate the proposed method.

## 1. Introduction

The use of multiscale models throughout engineering is nowadays ubiquitous. For example, fluid flow problems in heterogeneous media, the characterization of material properties such as conductivity, deformation or crack propagations in composite materials or chemical processes in biology all need mathematical models taking into account different physical processes at different scales. While in some applications different physical models might be used on different scales (quantum mechanics, molecular mechanics or continuum mechanics), we focus, in this contribution, on physical models described by partial differential equations (PDEs) with multiple scales.

Consider therefore a family of PDEs with appropriate boundary conditions
1.1
parametrized by *ε*, where and *Ω* is an open subset of . The parameter *ε* emphasizes the multiscale nature of the above family of PDEs, and represents a typical microscopic length scale of a heterogeneity in the system (oscillatory source term *f*^{ε} or multiple microscopic length scale considered by indexing the family of PDEs by *ε*=(*ε*_{1}(*ε*),…,*ε*_{N}(*ε*)) could be considered as well).

Numerous numerical techniques such as the finite-element method (FEM), the finite-difference method (FDM) or the finite-volume method (FVM) are nowadays available for the numerical discretization of (1.1). A major issue, however, is the need to resolve the finest length scale in the problem leading to a typical grid or mesh size *h*<*ε*. For example, thermal management of a multiphase composite with typical microstructures of about a hundred micrometres would need a mesh of several billion points for a three-dimensional computation with a material of size . Fortunately, appropriate averaging techniques such as the homogenization method have been developed in the past several decades [1–7] that allow for alternative numerical strategies. In these approaches, one characterizes the limit *u*^{ε} (or the limit of an appropriate flux of *u*^{ε}) as *ε* goes to zero and identifies an averaged equation *L*_{0}(*u*^{0})=*f* for the identified limit *u*^{0}. Several questions then arise such as whether the sequence *u*^{ε} converges and in which sense, whether the limit function (if it exists) solves a PDE, whether this PDE can be determined and finally whether *u*^{0} is a good approximation of *u*^{ε}. The rigorous treatment of these questions is at the core of the mathematical homogenization theory.

The homogenization theory is also at the core of most of the numerical methods for PDEs with multiple scales. We mention for example

— methods that supplement oscillatory functions to a coarse FE space, pioneered by Babuška & Osborn [8], generalized through the so-called multiscale finite-element method (MsFEM) [9], developed since then by many authors (MsFEM using harmonic coordinates [10,11], see [12] for a survey and additional references),

— methods based on the variational multiscale method (VMM) introduced in [13] and the residual free bubble (RFB) method [14] that are closely related to MsFEM-type strategy for homogenization problems [15],

— methods based on the two-scale convergence theory and its generalization [3,16] as proposed in [17] and developed in [18] using sparse tensor product FEM,

— projection-based numerical homogenization method based on projecting a fine-scale discretized problem into a low-dimensional space and eliminating successively the fine-scale components [19,20], and

— numerical homogenization methods that supplement effective data for coarse FE computation and approximate the fine-scale solution via reconstruction such as the heterogeneous multiscale method (HMM) [21,22] or related micro–macro methods [23–26].

In this paper, we focus on the aforementioned HMM. In the context of multiscale PDEs, this method relies on the following steps.

— A macroscale method such as the FEM, the FDM or the FVM defined on a macroscopic triangulation of the physical domain . The macroscale method solves an upscaled partial differential equation

*L*_{H}(*u*^{H})=*f*, where*L*_{H}is an*a priori*unknown approximation of*L*_{0}recovered from microscale computations.— Constrained microsimulations defined on a microscopic triangulation of sampling domains

*K*_{δ}_{j}=*x*_{K}_{j}+*δY*, where*Y*=(−1/2,1/2)^{d},*δ*≥*ε*, and*x*_{K}_{j}∈*K*are appropriate quadrature points. The microscale method solves a problem involving the original differential operator*L*_{ε}(⋅) usually with zero forces and with boundary conditions imposed from the macrostate*u*^{H}.

Our main aim is to present a reduced-order modelling technique, that can be combined with numerical homogenization techniques such as the HMM, to address the complexity issue of the classical numerical homogenization methods. This method, called the reduced basis (RB) finite-element heterogeneous multiscale method (RB-FE-HMM), was introduced in [27], combined with adaptive macroscopic methods in [28] and generalized for a class of nonlinear problem in [29]. In this contribution, we want to review the RB-FE-HMM and present a unified framework for this method by explaining its use for a variety of problems, including multiscale elliptic, parabolic and wave equations and a class of nonlinear elliptic or parabolic multiscale problems. In the RB-FE-HMM, a low-dimensional subspace, the so-called RB space, is constructed in an offline stage by a greedy algorithm. The offline stage is only performed once, and the outputs can be repeatedly used for many-query contexts in a so-called online stage, where the unknown effective parameters of the macroscopic solution are computed in the RB space. The moderate dimension of the RB space results in a computational cost for the online stage often comparable to a single-scale FEM. As demonstrated in the numerical examples, the RB-FE-HMM presents significant efficiency advantage over the FE-HMM and can be easily combined with different model problems and macrosolvers.

This paper is organized as follows. In §2, we introduce several model problems considered in this paper and briefly review the homogenization theory. The FE-HMM framework is reviewed in §3 for the various model problems and in §4 we present *a priori* error estimates and a complexity analysis for the FE-HMM. We present the RB-FE-HMM in §5 with a uniform description for all the model equations. The proposed reduced-order modelling strategy is then tested in §6 for several numerical examples, including wave propagation in inhomogeneous media and solute transport in unsaturated porous media.

## 2. Model problem, homogenization and finite-element heterogeneousmultiscale method

Here, we describe various PDEs with highly oscillatory coefficients and discuss briefly the averaging procedure called homogenization.

Our physical domain will always be a bounded polyhedron *Ω* in with *d*≤3. In order to explain our methodology, we consider second-order linear elliptic equations of the form
2.1
where *f*∈*L*^{2}(*Ω*). Here, we choose a zero Dirichlet boundary condition for the sake of simplicity.

We then explain how the same methodology can be applied to the following PDEs with appropriate initial and boundary conditions:

— linear parabolic equations 2.2

— linear wave equations 2.3

— nonlinear elliptic equations 2.4

— nonlinear parabolic equations 2.5

In equations (2.2)–(2.5), *a*^{ε} is a linear or nonlinear tensor that oscillates rapidly in space at the scale *ε*, which denotes a small scale in the problem such as the size of a typical heterogeneity under consideration in a porous medium, or the size of a typical microstructure in a composite material, etc. Solving any of the above equations by a standard FEM (or any other numerical method) requires, for small *ε*, a very fine mesh size *h*<*ε*, leading to a prohibitive computational cost.

### (a) Homogenization

In mathematical homogenization, one aims to describe an averaged equation corresponding to one of the class of PDEs with rapidly oscillating coefficients described previously. We describe briefly the homogenization procedure for the elliptic equation (2.1) and comment on similar techniques for the other equations. The formal approach based on asymptotic expansion consists of postulating an expansion *u*^{ε}(*x*)=*u*^{0}(*x*,*x*/*ε*)+*εu*^{1}(*x*,*x*/*ε*)+*ε*^{2}*u*^{2}(*x*,*x*/*ε*)+⋯ for the solution of (2.1). Here, we assume that *a*^{ε} is locally periodic, i.e. *a*^{ε}(*x*)=*a*(*x*,*x*/*ε*)=*a*(*x*,*y*) is *y*-periodic in *Y* (usually, *Y* is taken as the unit cube (−1/2,1/2)^{d}) and correspondingly, we assume that the functions *u*^{i}(*x*,*x*/*ε*)=*u*^{i}(*x*,*y*) are periodic in the second variable. Inserting the asymptotic equation in the original PDE and identifying the power leads to an averaged (homogenized) PDE depending on an averaged tensor *a*^{0}(*x*) that no longer depends on *ε*. For each macrolocation *x*, the explicit formulae of *a*^{0}(*x*) depending on the solution of a so-called cell problem (also referred to as a microproblem in this paper) are available for locally periodic problems. The solution *u*^{0}(*x*) of the homogenized PDE is called the homogenized solution [2]. To make this formal computation rigorous, one can use Tartar's method of oscillating test functions [30] (see also [2]) to show that , weakly in (*L*^{2}(*Ω*))^{d}.

Departing from the locally periodic case, there exists more general theory such as the *H*-convergence [30]. As a starting point, we have to consider a family of equations, corresponding to a family of tensors *a*^{ε} indexed by *ε*, that are uniformly elliptic and bounded, i.e. there exist positive such that for any
2.6
The *H*-convergence ensures then the existence of a subsequence of the matrices *a*^{ε} and a homogenized tensor *a*^{0} (again uniformly elliptic and bounded) such that for the corresponding subsequence, *u*^{ε} and *a*^{ε}∇*u*^{ε} converge weakly to *u*^{0} in and weakly to *a*^{0}∇*u*^{0} in (*L*^{2}(*Ω*))^{d}, respectively. For non-periodic oscillating tensors, the homogenized tensors *a*^{0}(*x*) are, in general, not known in an explicit form. Similar averaging procedure exists for the multiscale problem (2.2)–(2.5) (see [2,5,31,32]). For a practical solution, one has to rely on numerical approximation. The numerical methods that we describe below do not rely on periodic tensors. The scale separation seems required for our numerical strategy to make sense. However, the case of locally periodic tensors will sometimes be considered to derive a complete numerical analysis for the proposed multiscale method. Already in this situation, the determination of the homogenized tensor *a*^{0}(*x*) depends on the macrolocation *x*∈*Ω*, and we thus have an infinite number of cell problems to solve to obtain *a*^{0}(*x*) and a proper approximation of *a*^{0}(*x*) is required.

## 3. Numerical homogenization, micro–macro methods

We now present a numerical homogenization method, called the finite-element heterogeneous multiscale method (FE-HMM) [33,34], that is able to compute an approximation of the homogenized solution *u*^{0}(*x*), relying on a finite number of cell problems chosen in such a way that the overall computation is efficient and reliable. In a second step, an approximation of the fine-scale solution *u*^{ε} can be obtained by a reconstruction procedure.

### (a) Main ingredients

Assume for the sake of simplicity that *Ω* is a polyhedral domain in where *d*≤3 and consider a shape-regular family of partitions of *Ω* in simplicial or quadrilateral elements of diameter *H*_{K}, where we denote . As *H* is not required to be smaller than or even commensurate to *ε*, we call this triangulation a macroscopic triangulation of *Ω*. In its simplest form, the FE-HMM relies on the following ingredients:

(i) a macroscopic FE method based on a macroscopic triangulation of

*Ω*,(ii) a quadrature formula on each macroscopic element

*K*of the macroscopic triangulation, and(ii) microscopic FE methods defined on sampling domains around the integration points in

*K*used to recover the effective parameters (e.g. macroscopic conductivity) around the integration points.

We next describe the different ingredients listed above. A commonly used macroscopic FE space is given by
where is the space of polynomials on *K* of total degree at most ℓ if *K* is a simplicial FE, or the space of polynomials on *K* of degree at most ℓ in each variable if *K* is a quadrilateral FE. We note that for some problems for which mass conservation is required, other macroscopic FE space should be used. We mention, for example, the discontinuous Galerkin FE-HMM proposed and analysed for elliptic problem in [35] and for advection–diffusion problem (with possible high Peclet number) in [36].

For each element *K* of the macropartition, we consider a quadrature formula (*ω*_{Kj},*x*_{Kj})_{j=1,…,J} with weights *ω*_{Kj} and nodes *x*_{Kj} fulfilling classical assumptions (see §4 and [37]).

Finally, the microscopic FE method is defined as follows. Define for each quadrature node *x*_{Kj} a sampling domain *K*_{δj}=*x*_{Kj}+*δY* . On this sampling domain, we consider a simplicial micro-mesh and the micro-FE space *V* _{h}(*K*_{δj}) defined by
3.1
where the choice of *W*(*K*_{δj}) determines the boundary conditions used for computing the microfunctions . We consider two different spaces

— periodic coupling: and

— Dirichlet coupling: .

The microproblems depend on the specific problem (2.1)–(2.5). For (2.1), given a macroscopic function *v*^{H}∈*V* _{H}(*Ω*), we consider the linearization and the following problem: find such that and
3.2

### (b) The finite-element heterogeneous multiscale method for linear problem

At the macroscopic level, the numerical method is defined as follows: find *u*^{H}∈*V* _{H}(*Ω*) such that
3.3
where
3.4
In (3.4), (respectively ) denotes the solution of the microproblem (3.2). The following reformulation of the above macroproblem is useful for the analysis of the method, namely
3.5
where
3.6
with *τ*=(*x*_{Kj},*m*) and and is the solution of (3.2) that can be rewritten as
3.7
where **e**_{m}, *m*=1,…,*d* denotes the canonical basis of . For a proof of this equivalence, we note that can be represented by a linear combination of (see [38] for details). We next explain how the FE-HMM can be generalized for parabolic, wave and nonlinear equations. We focus only on the reformulation of the form (3.5). We close this section by noting that for practical computation *f*(*x*) should be replaced by an appropriate approximation *f*_{H}(*v*^{H}) in a FE space. Here and in what follows, we work with exact right-hand side for the sake of simplicity. The FE-HMM for parabolic homogenization problem (2.2) reads: find , such that
3.8
while for the wave equation it reads: find such that
3.9
where the bilinear form *B*_{H}(⋅,⋅) is defined by (3.5) for both problems and initial and boundary conditions must be supplemented.

### (c) The finite-element heterogeneous multiscale method for nonlinear problems

We start with the nonlinear homogenization problem (2.4). The FE-HMM reads: find *u*^{H}∈*V* _{H}(*Ω*) such that
3.10
where
3.11
The components of the tensor *a*^{0,h}(*x*_{Kj},*s*) are defined as in (3.6) using (3.7), where in both equations *a*^{ε}(*x*) must be replaced by *a*^{ε}(*x*,*s*). The parameter *τ* for now depends on *τ*=(*x*_{Kj},*s*,*m*).

Using Newton iterations for the nonlinear problem (3.10), we consider a sequence satisfying the following iteration scheme:
3.12
where ∂*B*_{H}(*z*^{H};*v*^{H},*w*^{H}):=*B*_{H}(*z*^{H};*v*^{H},*w*^{H})+*B*_{H}′(*z*^{H};*v*^{H},*w*^{H}), and
where ∂_{s}*a*^{0,h} denotes the derivative of *a*^{0,h} with respect to the second variable.

Finally, for nonlinear parabolic problems of the type (2.5) the FE-HMM reads: find , such that 3.13

## 4. *A priori* estimates, fine-scale reconstruction and complexity

Here, we discuss *a priori* error estimates for the various FE-HMMs introduced above. First, we observe that appropriate conditions on the quadrature formula are needed to ensure well-posedness of the macroproblem and optimal convergence. The following assumptions are the usual requirement for single-scale FEM with numerical quadrature (see [37]). We consider nodes and weights on the reference element (the corresponding nodes and weights on element are then obtained using a *C*^{1}- diffeomorphism from to *K*)

(Q1) , , ;

(Q2) , where if is a simplicial FE, or if is a rectangular FE.

### (a) *A priori* estimates

The general methodology to estimate the error between the FE-HMM solution *u*^{H} and the homogenized solution *u*^{0} of any of the problems (2.1)–(2.5) is based on a decomposition in macro-, modelling and microerrors as described below (see also [39]). We first introduce two auxiliary FE functions: *u*^{0,H}, the solution of any of the homogenized problems obtained by an FEM with numerical quadrature; and , the FE-HMM solution of any of the problems (2.1)–(2.5), but with a form obtained with exact microfunctions. Both *u*^{0,H} and are only introduced to analyse the various contribution to the error and we emphasize that these solutions cannot be obtained in practical applications. Indeed, to obtain *u*^{0,H}, one needs to know the exact homogenized tensor that is not available in general and to obtain one needs to know the exact solution of the microproblems. We then consider the following decomposition:
4.1
where ∥⋅∥ stands for the *L*^{2} or *H*^{1} norms for elliptic problems, the *L*^{2}([0,*T*];*H*^{1}(*Ω*)) or norms for parabolic problems, and the or norms for the wave problem. Of course, the rigorous analysis depends on the type of problem under consideration, but one common feature is that variational crimes are committed in the FE-HMM in the sense that the exact homogenized form *B*_{0}(⋅,⋅) differs from its numerical counterpart *B*_{H}(⋅,⋅) and hence standard Galerkin orthogonality arguments fail. This complicates the analysis especially for nonlinear problems [40,41], non-conforming FE discretization [35,36] or time-dependent problems [42–44].

Another common issue is that the microerrors are transmitted to the macroscale resulting in an error in the effective data. The so-called fully discrete analysis, first given in [33] for the FE-HMM, gives an indication of the complexity of the numerical method and indicates how to balance micro- and macro-mesh sizes in order to achieve a given accuracy with a minimal computational cost.

Finally, the modelling error encodes the geometric error owing to the mismatch between macrocomputational domain size and the size of the microperiod and the error incurred in imposing (artificial) boundary conditions in the microsampling domains (determined by the choice of the micro-FE space (3.1) that sets the coupling conditions between micro- and macro-FE functions). To be more specific, we mention the *a priori* error estimates for elliptic problems obtained in [33,34,39].

### Theorem 4.1

*Let u*^{0} *be the homogenized solution corresponding to the problem (2.1) and u*^{H} *be the solution of problem (3.3). Assume that (2.6), (Q1), and (Q2) hold. Then, under sufficient regularities of the tensor a*^{ε} *and the right-hand side f, we have the following estimates:*
4.2
*and*
4.3
*where C is independent of H,h and ε.*

The modelling error does not depend on the micro- or macro-mesh size. For locally periodic coefficients, it is possible to show that *e*_{mod}=0 using a slightly modified FE-HMM (collocated version) when and for the micro-FE space subset of . When and , one can show that *e*_{mod}≤*C*(*ε*/*δ*+*δ*) [34].

*A priori* error estimates similar to theorem 4.1 have been obtained for parabolic problems of the type (2.2) in [43], for nonlinear elliptic problems of the type (2.4) in [40,41] and for wave problems of type (2.3) in [42].

### (b) Fine-scale reconstruction

While the numerical methods described in §3*b*,*c* allow one to compute an approximation *u*^{H} of the homogenized solution *u*^{0}, they do not allow one to approximate the fine-scale solution *u*^{ε} of (2.1)–(2.5) in the energy norms. We note however that in the *L*^{2} norm *u*^{H} is still a good approximation of *u*^{ε} owing to the error estimate ∥*u*^{ε}−*u*^{0}∥≤*Cε* that holds for sufficiently regular problems [5].

Using the FE-HMM, we can nevertheless recover a fine-scale approximation of *u*^{ε} following a post-processing procedure inspired by Oden & Vemaganti [45]. Let *D*⊂*Ω* be the region of the computational domain where we want to approximate *u*^{ε} and assume that *u*^{H} has been computed on *Ω*.

We start by explaining the reconstruction procedure for the linear problem (2.1). Define then *D*⊂*D*_{η}⊂*Ω*, where dist(∂*D*,∂*D*_{η})=*η* and consider the following problem: find such that
The following error estimate is valid under appropriate regularity assumptions [34]:
For locally periodic homogenization problems, a simpler and cheaper reconstruction can be obtained. Indeed, consider, for example, the FE-HMM with piecewise linear macro-FE functions and the microfunctions *u*^{h}−*u*^{H} available in each *K*_{δ} (observe that for piecewise linear macro-FEM, there is only one sampling domain per macroelement). Next, extend these microfunctions periodically in *K* (we denote this extension by *u*_{h,K}) and consider the reconstruction
4.4
Assuming that and then the following error estimate holds [33,34]:
where *C* is independent of *H*,*h*,*ε* and denotes a broken semi-norm. This result is based on corrector results for homogenization problems, i.e. a function *u*^{1,ε} such that *u*^{ε}≃*u*^{0}+*u*^{1,ε}. For locally periodic problems, such function *u*^{1,ε} is obtained from the solution of localized microproblems of the type (3.2) (in a periodic Sobolev space), and the convergence results hold for smooth problems [5].

For the nonlinear problem (2.4), a similar procedure can be used [40,41], thanks to the results in [46], section 3.4.2. There it is shown that any corrector *u*^{1,ε} for the linear problem obtained from (2.4) by replacing *a*^{ε}(*x*,*u*^{ε}(*x*)) with *a*^{ε}(*x*,*u*^{0}(*x*)), where *u*^{0} is the solution of the corresponding homogenization problem, is also a corrector for the solution *u*_{ε} of the nonlinear problem (2.4) and that
4.5
Hence, a similar reconstruction as defined in (4.4) can be used and for locally periodic problems assuming again and we have
where *C* is independent of *H*,*h*,*ε* and *r*_{ε} is defined in (4.5).

### (c) Complexity

The convergence rates in theorem 4.1 show a classical rate for the macroerror and a better than usual rate in the microerror (loosely speaking, this is due to the fact that the product of microfunctions enters in the bilinear form *B*_{H}(⋅,⋅); see [33,41,47] for details).

We next discuss the computational cost of the FE-HMM. If we denote *N*_{mic} the number of degrees of freedom (DOF) in each space dimension for the discretization of the sampling domain *K*_{δj}, we obtain *h*=*δ*/*N*_{mic} and hence . By noting that *δ* scales with *ε* (e.g. *δ*=*Cε* with *C* a constant of moderate size), we have . We next denote by the number of DOF for the micro-FEM and by *M*_{mac} the number of DOF of the macro-FEM. The macro-mesh size *H* and the micro-mesh size are related to *M*_{mac} and *M*_{mic} (for quasi-uniform macro-meshes) as
and according to the *a priori* error estimates of theorem 4.1 optimal macroscopic convergence rates require
with an induced complexity given by
where *n*_{s} denotes the number of sampling domains per macroelement .

We first note that the method is *independent* of the size of the oscillatory parameter *ε*. This is in sharp contrast to classical FEMs that would require a computational cost that scales with *ε* as . Second, it can be seen, as first noted in [33], that the complexity is superlinear with respect to the macro-DOF. As an example, if we choose piecewise linear simplicial FEs and assume that the complexity is proportional to the total DOF, then we obtain a cost of (*H*^{1} norm) and (*L*^{2} norm).

Of course the method is well suited for parallel implementation as the microproblems are solved independently. We also note that the memory requirement is proportional to *M*_{mac}+*M*_{mic} only as the microproblems, being independent of one another, can be solved one at a time. Finally, as investigated in [48], using spectral method or p-FEM for the microsolvers can reduce the complexity of the FE-HMM down to a log-linear complexity. This approach however requires high regularity of the oscillating tensor *a*^{ε}. Another tool allowing for a reduction of the computational cost is the use of adaptive techniques at the macroscale. Indeed, microcomputations can be recycled in macroelements that are not refined [49].

## 5. Reduced-order modelling numerical homogenization

The key issue that leads to a superlinear computational cost in a numerical homogenization method such as the FE-HMM is the need of repeated computation of microproblems (3.7) at each quadrature point. At the same time in view of theorem 4.1, an increasing number of micro-DOF as the macroscopic mesh gets refined is needed. This issue has triggered the development of a reduced-order modelling strategy for the FE-HMM. The framework is built on the so-called RB methodology [50–52] first used in the context of numerical homogenization in [53,54] and for the FE-HMM in [27–29]. The new method is called the RB finite-element heterogeneous multiscale method (RB-FE-HMM).

The main observation is the following: instead of computing microproblems in each macroelement at the quadrature points, we identify in *an offline stage* a small number *N* of pre-computed representative microsolutions to construct a RB space. The selected parameters *τ* which determine those representative microsolutions are selected by a greedy algorithm based on a large parameter training set (figure 1). A key tool in the greedy algorithm is the use of an *a posteriori* error estimator to select the parameter *τ* for which the microsolutions vary the most.

The actual macrocomputation is performed in an online stage, and the missing effective data are computed at the required quadrature points of the macroelements by solving the microproblem in the pre-computed RB space (of dimension *N*) which leads to solving small linear systems of size *N*×*N*.

We see that for the RB-FE-HMM, the repeated micro-FEM computation and micro-mesh refinement are avoided though the pre-computation of a fixed low dimensional approximation space (the RB space). In addition, the RB space is independent of the macrosolvers which can be repeatedly used with different macroscopic partitions and different types of problems (static or evolutionary). In what follows, we give some details of both the online and the offline stages.

### (a) Offline stage

We observe that the solutions of the microproblems for any of the problems (3.3), (3.9), (3.10) or (3.13) depend on the parameter *τ*=(*x*,*m*) (linear problems) or *τ*=(*x*,*m*,*s*) (nonlinear problems). To treat both cases at the same time, we set *τ*=(*κ*,*m*), where *κ*=*x* (linear case) and *κ*=(*x*,*s*) (nonlinear case). The fundamental condition for the efficiency of the RB method is that the multiscale tensor *a*^{ε}(*x*,*s*) has an affine representation. We first set a correspondence of an arbitrary sampling domain included in *Ω*, namely *K*_{δ}=*x*+*δy* with *y*∈*Y* =(−1/2,1/2)^{d} through the affine transformation
5.1
We can then map the tensors *a*^{ε}(*x*) or *a*^{ε}(*x*,*s*) into the reference domain *Y* =(−1/2,1/2)^{d}
5.2
When it yields no confusion, we simply write *a*_{κ}(*y*) for either situation. A crucial assumption for the RB methodology is that *a*_{κ}(*y*) has an affine representation of the form
5.3
For example, for a tensor of the form , where *I* is the *d*×*d* identity matrix and *α* is a scalar, the above representation exists. When such an explicit representation is not available, one can use a greedy algorithm, called the empirical interpolation method, to approximate a non-affine tensor by an affine one of the form (5.3) (see [55]).

We next map the microproblem (3.7) (or its nonlinear version) for an arbitrary sampling domain *K*_{δ}=*x*+*δy* included in *Ω* to the reference domain *Y* by using the change of variable *G*_{x}(*y*)
5.4
where *τ*=(*κ*,*m*) and .

We next choose a set of parameters where is a compact subspace of *Ω* for linear problems and a compact subspace of for nonlinear problems. The construction of the RB space is carried out in an offline procedure, where a small number {*τ*_{1},…,*τ*_{N}} of representative microproblems with are selected by a greedy algorithm. The corresponding FE solutions of the cell problem (5.4) for these selected parameters will give the RB functions obtained from through normalization and orthogonalization. Hence, we obtain the RB space .

### Algorithm 5.1

*(Offline procedure.) Define a training set Ξ*_{RB}, *an offline tolerance tol*_{RB}, *and compute for a random selected τ*_{1}∈*Ξ*_{RB} *the first RB function*

(

*i*)*assume that**is computed, compute**where**is the microsolution solved in S*_{ℓ};(

*ii*) if*the offline stage ends (appropriate data must then be stored, see (5.6) and (5.8)); otherwise, continue with the next step*;(

*iii*)*select the next representative parameter**and compute the new basis function*.*Let ℓ=ℓ+1 and go back to step*(*i*).

The offline procedure is operated only once, and the outputs can be repeatedly used for the later online computation. Therefore, we require the micro-FEM used in the offline stage to be very accurate, so that the corresponding micro-FE error *e*_{mic} does not affect the online results.

We note that the direct computation of in algorithm 5.1 can become quite expensive owing to the computation of highly resolved micro-FE solution over a large training set. This issue can be resolved by estimating by an *a posteriori* error estimator Δ_{τ} which can be computed by solving a few pseudo-FE solutions. The *a posteriori* estimator for linear problems is designed in [27]. For nonlinear problems, some care is needed to construct such estimators which have been derived in [29] in a numerical homogenization context. The following results have been obtained:
We also emphasize that the *a posteriori* error can be used in the online stage to certify the accuracy of the online solution. Appropriate procedures to compute the constant *C* are available. For nonlinear problems, Newton iteration is applied and we thus need to have control on the derivative of with respect to the parameter *s*. A general result proved in [56,57] shows that if the best *N*-dimensional approximation of a subset of a Hilbert space has a rapidly decaying projection error (e.g. exponential decay), then the *N*-dimensional space obtained from the RB method enjoys the same projection error.

### (b) Online stage

Compared with the FE-HMM, the advantage of the RB-FE-HMM is that in the online stage, all the microproblems are solved in the same RB space *S*_{N} with dimension *N* (usually of moderate size). Thus, the computational cost for the online stage is and scales *linearly* with the macroscopic DOF (assuming that the cost is proportional to the total DOF).

We next detail the online procedure. To obtain the numerical homogenized solution, we still need to solve the macroproblem either in form of (3.3) or (3.10). Now, the unknown homogenized tensor *a*^{0}(*κ*) can be estimated by
5.5
where . By expending and using the affine representation (5.3), we can further write (5.5) as
where *β*_{τ}:=(*β*_{1,τ},…,*β*_{N,τ}). The matrices *F*_{p,n} and *G*_{p} are the offline outputs defined as
5.6
for *m*,*n*=1,…,*d*.

In order to obtain the coefficients *β*_{τ}, we need to solve the following cell problem:
5.7

Using (5.6), (5.3) and the following offline output matrices for *p*=1,…,*P*,
5.8
equation (5.7) can be written as an *N*×*N* linear system
5.9
We note that the necessary data that need to be stored at the end of the offline stage are only a few *N*×*N* matrices with low storage requirement. The linear system (5.9) is independent of the macropartition or the macrosolvers. Furthermore, for linear multiscale problems (5.9), it is independent of the right-hand-side function *f* in the multiscale model equations nor does it depend on the boundary and initial conditions for time-dependent problems. Therefore, the offline outputs can be used for various macroscale scenarios.

In comparison, the microcell problems for the FE-HMM (3.7) are solved by the micro-FEM with the number of DOF *M*_{mic}≈*M*_{mac}. This triggers simultaneous micro- and macrorefinement and leads to a significant computational overhead compared with the RB-FE-HMM. We study this significant difference in the performance of the FE-HMM versus RB-FE-HMM in the numerical example of §6.

Finally, we end this section by mentioning that the *a priori* error estimate of the RB-FE-HMM is similar to the FE-HMM but with one more error term arising from the RB model reduction, i.e.
5.10
where *u*^{H,RB} is the solution of (3.3) or (3.10) based on the RB-FE-HMM. Compared with the FE-HMM, the errors *e*_{mac} and *e*_{mod} remain the same, but *e*_{mic} is much smaller owing to the offline requirement that the representative microproblems are computed accurately. The term *e*_{RB} is bounded by the *a posteriori* estimator and therefore controlled by the given offline tolerance. This *a posteriori* error estimator can also be computed in the online stage to certify the accuracy of the computed macrosolution by quantifying the RB error. *A priori* estimate for *e*_{RB} relies on appropriate assumption on best *N*-dimensional subspace that minimizes the projection error of arbitrary functions in the space of solutions of cell problems (see [27,29]). In practical computations, the RB-FE-HMM and the FE-HMM have nearly the same accuracy, provided *simultaneous refinement* of macro- and micro-meshes is implemented for the FE-HMM. By contrast, the computational cost for the online stage of the RB-FE-HMM is comparable to the cost of a FEM for single-scale problems and yields a substantial saving when compared with numerical homogenization methods such as the FE-HMM.

## 6. Numerical examples

Here, we illustrate the performance of the RB-FE-HMM for three numerical examples related to various models (2.1), (2.3), (2.5) presented in §2. In the first example, we consider a stationary linear elliptic multiscale problem with a tensor displaying discontinuity on the microsampling domains. The simulation of wave propagation in an inhomogeneous macrodomain is treated in the second problem. In the last example, we consider a Richards equation (a nonlinear problem) in a three-dimensional heterogeneous medium, a widely used model to evaluate the pressure head in soil infiltration models.

*Computational settings.* In the following numerical tests, we use a simplicial partition of the computational domain and piecewise linear polynomial basis functions. The quadrature points for the corresponding macro-FEM are at the bary centres of the elements. All the tests are performed in a single thread Matlab environment, based on the code presented in [58]. The numerical convergence rates are presented in relative errors, for example, ∥*u*^{H}−*u*^{0}∥/∥*u*^{0}∥.

### (a) Two-dimensional stationary problem with discontinuity on the microdomain

We consider in the macrodomain *Ω*=[0,1]^{2} a stationary problem of the form
with a diagonal locally periodic multiscale tensor with discontinuities illustrated in figure 2. We set *f*(*x*)=1 and pose a mixed boundary condition, i.e. *u*^{ε}(*x*)=0, *x*∈{*x*_{1}=0}∪{*x*_{1}=1}, the normal derivative ∂*u*^{ε}/∂*n*=0, *x*∈{*x*_{2}=0}∪{*x*_{2}=1}.

*Offline stage.* We conduct the offline process following algorithm 5.1. The offline settings and outputs are shown in table 1. It can be seen that only 10 RB functions are needed to reach our given precision tolerance of 10^{−10}. At the end of the offline stage, necessary data (5.6), (5.8) are stored for later online use, as discussed in §5*b*. Owing to the discontinuity of the affine representation in this example, we need to apply the so-called successive constraint method [51] to estimate the coercive factor *α*.

*Online stage.* Next, we apply the RB-FE-HMM online procedure to obtain the numerical homogenized solution *u*^{H,RB}. Here, we use uniform macro-meshes with sizes 16×16,32×32, 64×64, 128×128 and 256×256, respectively. The values of the unknown homogenized tensor on the macroquadrature points are estimated using the RB obtained from the offline stage for any macrodomain partition. In this test, we compare the RB-FE-HMM solution *u*^{H,RB} with the FE-HMM solution as the reference solution for *u*^{0} computed with a uniform 512×512 mesh as both macro- and micro-meshes. In figure 3, we show a log–log plot of the error ∥*u*^{H,RB}−*u*^{0}∥ in the *H*^{1} and *L*^{2} norms versus where *M*_{mac} is the macroscopic DOF. According to the RB-FE-HMM *a priori* error analysis [27], the errors in the *H*^{1} and *L*^{2} norms decay with macrorates and as confirmed by the numerical decay rates shown in figure 3. As discussed in §5*b*, the microerror *e*_{mic} of the RB-FE-HMM is which is in our setting and remains unchanged during all the online procedures, whereas the RB error *e*_{RB} is bounded by the given tolerance tol_{RB}. Therefore, both *e*_{mic} and *e*_{RB} are negligible compared with *e*_{mac}.

In table 2, we present the CPU time comparison between the RB-FE-HMM and the FE-HMM. As we can see, when the macro-mesh size reaches 256×256, the RB-FE-HMM online time is only 3% of the FE-HMM cost and even considering the offline overhead, the RB-FE-HMM is still more efficient than the FE-HMM. We conclude from this test that the RB-FE-HMM proceeds with a significant computational speedup compared with the FE-HMM especially for a relatively fine macro-mesh.

### (b) Wave propagation in inhomogeneous media

In this example, we report numerical performance of the RB-FE-HMM for the linear multiscale wave equation (2.3). Here, we consider a multiscale tensor that displays different heterogeneity in three subdomains as shown in figure 4*a*, where we marked each subdomain with one colour and the multiscale tensor is diagonal and set as follows:
We note in figure 4*b* that we have sharp media discontinuities across the different subdomains. We consider *u*^{ε}(*x*,0)=0.1*e*^{−((x1−1)2+(x2−1)2)/σ2}, with *σ*=0.1 as the initial condition and use zero Dirichlet boundary condition for the sake of simplicity.

*Offline stage.* Owing to the different media, we perform the offline procedure in each subdomain and obtain three sets of offline basis functions which will be used in the online stage to estimate the unknown data in the corresponding subdomains. The offline settings and outputs for this test are presented in table 3.

*Online stage.* Using the offline outputs in the online procedure, we obtain the RB-FE-HMM solution *u*^{RB} shown in figure 5*a,c,e* at time *t*=0,0.1, and 0.2, respectively. For the initial condition at time *t*=0, we choose a Gaussian pulse. When the wavefront reaches the subdomain *Ω*_{3}, a significant scattering occurs owing to the small homogenized diffusion coefficient in *Ω*_{3}. We next compare our multiscale algorithm with a numerical solution by using the local arithmetic average of the multiscale tensor *a*^{ε}(*x*) in each subdomain (figure 5*b,d,f*). We can observe that the profiles of the two solutions differ significantly.

We finally compare the RB solution *u*^{H,RB} to the standard FE solution obtained by solving the homogenized equation with the corresponding explicit homogenized tensor. The errors and are shown in figure 6. We can observe that the errors decay with rates and , which corroborates the analysis in [42].

### (c) Richards equation in an unsaturated soil domain

In the last example, we study a nonlinear parabolic multiscale problem similar to (2.5), known as the Richards equation in subsurface flow modelling. The classical Richards equation models the flow pressure head in unsaturated media and is often combined with a mass transportation equation in order to simulate pollutant distribution in unsaturated soil [59]. In the literature on numerical simulations for the Richards equation, the main study and experiments are usually carried out for single-scale problems. However, the soil often has a multiscale structure owing to the large scale range between the macrocomputational domain and the small pore structure in the soil. The solution of a three-dimensional multiscale Richards equation is thus a challenging task, and we next discuss computational results obtained with the RB-FE-HMM.

We consider the following Richards equation in a three-dimensional computational domain as presented in figure 7:
where *u*^{ε} is the pressure head of the pollutant flow, and the water content function *Θ*(*u*) is defined as
where the residual water content *Θ*_{r}, saturated water content *Θ*_{s} and *α*,*n*,*m* are model parameters shown in table 4.

We next assume that the multiscale tensor is diagonal with entries defined as where

In order to show the variation of the pressure more clearly, we define the relative pressure head as *v*^{ε}(*x*)=*u*^{ε}−*x*_{3} and therefore *v*^{ε}(*x*) satisfies
6.1
We assume that the initial condition is *v*^{ε}(*x*,0)=0.1*x*_{3}−0.5 and that on the top boundary ∂*Ω*_{top} of the domain *Ω* *v*^{ε} satisfies Dirichlet boundary condition *v*^{ε}(*x*,*t*)=0.1*x*_{3}−0.5, *x*∈∂*Ω*_{top}. Homogeneous Neumann boundary conditions are set on all the other boundaries. Homogenization for such problems was studied in [61], where it is shown that the homogenized Richards equation has the same form as (6.1) with an oscillating tensor replaced by a homogenized one.

*Offline stage.* As mentioned in §5*a*, for this type of nonlinear problem, we have cell problems defined in (5.2) indexed by parameter *τ*=*κ*,*m*, where *κ*=(*x*,*s*). We apply the RB-FE-HMM to (5.4) with *κ*=(*x*,*s*). The offline settings and outputs can be seen in table 5, where the parameter range of *s* is estimated by an offline procedure proposed in [62].

*Online stage.* The total evolution time we consider in this experiment is *t*=3 and we set the time step Δ*t*=0.005 for the time integrator. We use the linearized Picard scheme proposed in [63] and a macro-mesh with 12 304 DOF. The total RB-FE-HMM online computational time is 4197 s for 600 time steps and for each time step the CPU time cost is about 6.5 s. As for the CPU time comparison, we performed a computation with single-scale FEM with the same macro-mesh, for which each time step computation took around 2 s. Therefore, we can conclude that the RB-FE-HMM online time cost is comparable to a single-scale FEM (up to some constant). In figure 8, we show the evolution of *v*^{0}(*x*,*t*) from *t*=0 to the final time *t*=3, illustrating the variation of the pressure head in the soil domain.

## 7. Conclusion

We have presented a unified framework to combine reduced-order modelling techniques with numerical homogenization methods for a variety of problems, including linear and nonlinear elliptic, parabolic and wave equations with highly oscillatory data. The reduced-order modelling technique, built on the RB method, allows one to pre-compute a representative number of microfunctions in an offline stage. This pre-computed RB is then used in an online stage to compute effective data for a homogenized model at arbitrary locations in the computational domain. The use of RB also allows the outputs of the offline stage to be repeatedly used for many-query contexts. The accuracy of this representative basis is controlled by appropriate *a posteriori* error estimators in the offline stage that also permit one to certify the accuracy of the online solution. We have shown that two issues in the numerical approximation of PDEs with multiple scales are addressed by the RB-FE-HMM, namely

— a computational cost

*independent*of the size of the oscillatory parameter*ε*, thanks to the numerical homogenization techniques and— a linear computational cost for the online stage obtained through the RB available for the online microproblems.

Finally, we have tested the numerical method on a variety of problems, including an elliptic problem with discontinuous microscopic oscillatory data, wave propagation in inhomogeneous media, and a three-dimensional infiltration problem in an unsaturated porous medium. Substantial computational savings are observed for the new reduced-order modelling numerical homogenization method when compared with classical numerical homogenization.

## Acknowledgements

The research of A.A. and Y.B. is partially supported by the Swiss National Foundation 200021_134716.

## Footnotes

One contribution of 13 to a Theme Issue ‘Multi-scale systems in fluids and soft matter: approaches, numerics and applications’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.