Reduced basis ﬁnite element heterogeneous multiscale method for quasilinear elliptic homogenization problems

A reduced basis ﬁnite element heterogeneous multiscale method (RB-FE-HMM) for a class of nonlinear homogenization elliptic problems of nonmonotone type is introduced. In this approach, the solutions of the micro problems needed to estimate the macroscopic data of the homogenized problem are selected by a Greedy algorithm and computed in an oﬄine stage. It is shown that the use of reduced basis (RB) for nonlinear numerical homogenization reduces considerably the computational cost of the ﬁnite element heterogeneous multiscale method (FE-HMM). As the precomputed microscopic functions depend nonlinearly on the macroscopic solution, we introduce a new a posteriori error estimator for the Greedy algorithm that guarantees the convergence of the online Newton method. A priori error estimates and uniqueness of the numerical solution are also established. Numerical experiments illustrate the eﬃciency of the proposed method.


Introduction
Quasilinear elliptic problems enter the modeling of numerous problems such as phase transitions, flow in porous media, or reaction and diffusion in electrolysis to mention a few examples [12]. Numerical approximations of such problems have been analyzed by many authors. We mention the works of Douglas and Dupont [28], and Nitsche [32], where the first a priori error analysis was given for the finite element method (FEM). Much recently and relevant for the present work, we mention the analysis obtained in [11] for a FEM with numerical quadrature, i.e., when the continuous variational form originating from the nonlinear problem is approximated by a quadrature formula. In this paper we are interested in quasilinear elliptic problems with highly oscillatory data of the form −∇ · (a ε (x, u ε (x))∇u ε (x)) = f (x) in Ω, (1) in a domain Ω ⊂ R d , d ≤ 3, where a ε (x, u) = (a ε mn (u, s)) 1≤m,n≤d is a d × d tensor, associated to ε > 0, a sequence of positive real numbers going to zero and f ∈ H −1 (Ω). For simplicity we assume homogeneous Dirichlet boundary conditions u ε = 0 on ∂Ω but we emphasize that more general boundary conditions could be considered. Such problems arise for example in infiltration of water in an unsaturated porous media modeled by the (stationary) Richards equation [15] or (stationary) heat conduction in a composite material [30]. For the standard FEM, a finescale resolution is needed for a satisfactory approximation. If the ratio between the scale of interest and the finest scale in the problem is too large, the FEM approximation will have a prohibitive number of degrees of freedom (DOF), leading to an enormous computational cost. For efficient numerical computations, an appropriate upscaling of equation (1) is thus needed. Such coarse graining procedures are rigorously described by the mathematical homogenization theory [16,29] and are studied for the class of problems (1) in [13,18,27]. These analyses show that the solution u ε of (1) converges in a weak sense to u 0 as ε → 0, where the homogenized function u 0 is the solution of an effective (homogenized) equation that is of the same quasilinear type as the original equation with an effective homogenized tensor a 0 (x, u 0 (x)) that depends nonlinearly on u 0 . Numerical homogenization methods for problems of the type (1) are derived in [22] for the multiscale finite element method (MsFEM) and in [26,9] for the finite element heterogeneous multiscale method (FE-HMM) [5,25]. The MsFEM is based on a standard FE space enriched with oscillatory functions, while the FE-HMM is based on a strategy first proposed in [25] that consists in macroscopic FEM on a macroscopic mesh with quadrature formula (QF), with effective data (the homogenized tensor at the quadrature points) recovered on the fly from micro problems. These micro problems, defined on sampling domains centered at the macroscopic quadrature points of the QF, use only the oscillatory data given by the problem (1). We focus on the FE-HMM proposed in [26,9] for quasilinear problems. The practical implementation relies on a Newton method for the macroscopic nonlinear FEM. Since the value of the corresponding macroscopic solution is updated at each Newton iteration, the microscopic problems in each element of the macroscopic mesh need to be recomputed. Although the micro problems can be solved independently in parallel, the cost of the procedure mentioned above can be prohibitive, especially for high dimensional problems.
In this paper, we show how the use of the reduced basis (RB) method (see [34,33,35] and references therein) for computing the micro problems permits to considerably improve the efficiency of the standard nonlinear FE-HMM. The use of RB for numerical homogenization problems has first been proposed by in [19] and analyzed for the FE-HMM for a class of linear elliptic problems in [3]. The algorithm proposed in this paper for nonlinear problems relies on online and offline procedures: in the offline procedure accurate micro solutions for the original problem on sampling domains are selected and computed. These micro problems are parametrized by the location of the cell problem in the domain Ω and the macroscopic solution at this location. A greedy algorithm allows to choose an optimal basis of micro functions (computed with high accuracy) for selected values of the parameters. In the online stage, a Newton method for the RB-FE-HMM implementation is proposed with microscopic solutions computed in the reduced basis space, which amounts to solve small dimensional linear systems in each element of the macroscopic mesh. The overall computational cost of the online macroscopic Newton method is similar to the cost of single scale nonlinear problems. One difficulty is the design of an a posteriori error estimator in the offline stage that is both efficient and also guarantees that the online Newton method converges. We propose in this paper such a posteriori error estimators and prove the convergence of the online Newton method and the uniqueness of the numerical solution. Furthermore, a fully discrete error analysis of the quasilinear RB-FE-HMM is derived. This paper is organized as follows. In Sect. 2, we briefly recall the framework of homogenization theory in our context of quasilinear elliptic problems of nonmonotone type. We then present in Sect. 3 the new nonlinear RB-FE-HMM with its offline and online procedures, and analyze its convergence in Sect. 4. We explain some implementation issues in Sect. 5. Finally, numerical experiments in Sect. 6 show how the use of reduced basis considerably improves the efficiency by reducing drastically the number of degrees of freedom for various problems.

Homogenization of quasilinear elliptic problems
We assume that the tensor a ε (x, s) in (1) is uniformly elliptic and bounded with respect to s and ε, i.e., there exist λ, Λ 1 > 0 such that and that the functions a ε mn (x, s), m, n = 1, . . . , d are continuous, bounded and uniformly Lipschitz continuous with respect to s. Then, for all fixed ε > 0, the weak form of (1) has a unique solution u ε ∈ H 1 0 (Ω) (we refer for example to [23,Theorem 11.6] for a proof). The solution, for each ε, satisfies the a priori bound u ε H 1 (Ω) ≤ C f H −1 (Ω) , hence one can apply standard compactness arguments to the sequence of solution u ε that ensure the existence of a subsequence of {u ε } converging weakly in H 1 (Ω). The homogenization result is shown in [18,Theorem 3.6] (see also [27]) and reads as follows: there exists a subsequence of {a ε (·, s)} (again indexed by ε) such that the corresponding sequence of solutions {u ε } converges weakly to u 0 in H 1 (Ω). The limit function u 0 is the solution of the homogenized problem The tensor a 0 (x, s), called the homogenized tensor, can be shown to be Lipschitz continuous with respect to s, uniformly elliptic, and bounded [18,Prop. 3.5], i.e., there exists Λ 2 > 0 such that 1 a 0 (x, and there exist λ, Λ 1 > 0 such that a 0 satisfies (2) (possibly with different constants). Under these assumptions, the homogenized problem (3) has also a unique solution u 0 ∈ H 1 0 (Ω). We mention that for a locally periodic tensor of the form a ε (x, s) = a(x, x/ε, s) where a(x, y, s) is Y periodic with respect to y, the weak convergence of u ε to the solution of (3) holds for the whole sequence {u ε } and the homogenized tensor can be characterized in the following way [13]: where J χ(x,y,s) is a d×d matrix with entries J χ(x,y,s) ij = (∂χ i )/(∂y j ) and χ i (x, ·, s), i = 1, . . . , d are the unique solutions in W 1 per (Y ) := {z ∈ H 1 per (Y ); Y zdx = 0} of the linear cell problems with parameters x ∈ Ω, s ∈ R Y a(x, y, s)∇ y χ i (x, y, s) · ∇w(y)dy = − Y a(x, y, s)e i · ∇w(y)dy, ∀w ∈ W 1 per (Y ), (6) where H 1 per (Y ) := {g ∈ H 1 (Y )| g periodic in Y } and e i , i = 1, . . . , d are the vectors of the canonical basis of R d .
Remark 2.1. We sometimes refer to the problems (3) or (1) as "non monotone problems". This stems from the following fact: writing for example (3) in weak form we observe that the monotonicity property with C ≥ 0 does not hold in general for the quasilinear problem (3) (or (1)). This lack of monotonicity makes the numerical analysis for FEM a nontrivial task, in particular when quadrature formula are used [11].
For our analysis, we will further assume that the tensor a ε is symmetric (and thus also a 0 ) and that the homogenized tensor is continuous, 3 Reduced basis FE-HMM for quasilinear problems As the homogenized tensor a 0 in (3) is in general unknown, the task in numerical homogenization is to design an algorithm capable of computing an approximation of the homogenized solution u 0 without knowing a 0 , relying on a finite number of localized micro problems, i.e. cell problems, chosen in such a way that the overall computation is both efficient and reliable.
Here, we generalize the RB-FE-HMM introduced in [4] for linear elliptic problems to quasilinear elliptic problems. This method relies on a macroscopic solver with macroscopic data recovered by microscopic simulations (the micro problems) performed on sampling domains located at appropriate quadrature points of the macroscopic mesh. In addition, in order to avoid repeated micro computations, the solution of the micro problem are computed in finite dimensional space of low dimension spanned by a so-called reduced basis obtained in an offline procedure.

Preliminaries
We describe here the macro and micro finite element spaces needed to define and analyze the RB-FE-HMM.
Macroscopic mesh and quadrature formulas. The RB-FE-HMM is based on a macro finite element (FE) space where T H is a shape-regular family of (macro) partition of Ω in simplicial or quadrilateral elements K of diameter H K , and R (K) is the space P (K) of polynomials on K of total degree at most if K is a simplicial FE, or the space Q (K) of polynomials on K of degree at most in each variable if K is a parallelogram FE. For a given macro partition, we define as usual H := max K∈T H H K . We highlight that H in our discretization is allowed to be much larger than ε.
For each element K of the macro partition we consider an affine transformation F K such that K = F K (K), whereK is the reference element (simplicial or parallelogram). For a given quadrature formula {x j ,ω j } J j=1 onK, the quadrature weights and integration points on K ∈ T H are then given by ω K j =ω j |det(∂F K )|, x K j = F K (x j ), j = 1, . . . , J. We make the following assumptions on the quadrature formulas, which are standard assumptions also for linear elliptic problems [24]: Microscopic mesh and RB. We consider a micro FE space S q (Y, Tĥ) ⊂ W (Y ) with simplicial or quadrilateral FEs and piecewise polynomial of degree q on the domain Y = (−1/2, 1/2) d equipped with a conformal and shape regular family of triangulation denoted Tĥ. The space W (Y ) denotes either the Sobolev space for a periodic coupling or for a coupling with Dirichlet boundary conditions. We then consider the RB space, which is a subspace of S q (Y, Tĥ) with a low dimension denoted S N (Y ) = span{ξ n,N (y), n = 1, .., N } ⊂ S q (Y, Tĥ).
whereξ n,N (y), n = 1, .., N denotes the reduced basis. Notice that for the analysis of the RB-FE-HMM, we shall also consider a RB space of the form S N (Y ) = span{(ξ n,N ,ζ n,N ), n = 1, .., which is a subspace of dimension N of (S q (Y, Tĥ)) 2 involving the same functionsξ n,N as in S N (Y ) and whereζ n,N ∈ S q (Y, Tĥ), n = 1, .., N . The construction of the RB spaces S N (Y ) and S N (Y ) is discussed in Sect. 3.4 below. For each macro element K ∈ T H and each quadrature point x K j ∈ K, j = 1, . . . , J, we define the sampling domains K δ j = x K j + (−δ/2, δ/2) d , (δ ≥ ε). We observe that each sampling domain K δ j is in correspondence with Y through the affine transformation This transformation applied to the RB space (10) permits to define the RB space S N (K δ j ) associated to each sampling domain K δ j as

Online procedure: the RB-FE-HMM
Assuming that the RB space has been pre-constructed in the offline stage described in the next section, we introduce a macro method similar to the FE-HMM with the micro problems solved in the RB space. The nonlinear RB-FE-HMM for (1) is defined as follows: find u H,RB ∈ S 0 (Ω, with a bilinear form defined for all u H , v H , w H ∈ S 0 (Ω, T H ) by where for the scalar parameter and similarly for w s N,K j (x). The problem (15) requires the solution of an N ×N linear system, where the details of the offline output and the online implementation are discussed in Sect. 5. The efficiency of the RB procedure relies in the fact that the dimension N of the RB space is usually small. Furthermore, in contrast to the standard FE-HMM, the number of degrees of freedom (DOF) of the micro (RB) space remains fixed during the online procedure and does not increase as the macroscopic DOF increase. This is in sharp contrast with the FE-HMM for which the simultaneous refinement of the macro and micro DOF is a major computational issue [1].

Solution of the macro quasilinear problem and Newton method
While the cell problems (15) are linear, the macroscopic problem (14) is nonlinear and is usually solved by a Newton method.
The following reformulation of the bilinear form of the RB-FE-HMM will be useful to define the Newton method used in practice to compute a numerical solution u H,RB of (13). The bilinear form (14) can be rewritten as where we define the numerical homogenized tensor as . . , d is the solution of a cell problem (see (28) below) corresponding to the sampling domain K δ j .
Inspired by [28,9], we explain here how to solve the nonlinear problem (13) with the Newton method. For given z H , v H , w H ∈ S 0 (Ω, T H ) we first define the Fréchet derivative ∂B H obtained by differentiating the nonlinear quantity where by the reformulation of the RB-FE-HMM bilinear form (16) we derive The Newton method for approximating a solution u H of the nonlinear RB-FE-HMM (13) by a sequence {u H k } reads in weak form The fact that the Newton method is well defined and convergence is discussed in Sect. 4.2 while an efficient implementation is detailed in Sect. 5.

Offline procedure: RB for quasilinear problems
This section describes the offline stage of the RB algorithm in our context of quasilinear elliptic problems. The task is to construct a low dimensional RB space S N (Y ) spanned by a small number N N of representative solutions of the cell problems (25) below (depending on the quadrature node s K j and the nonlinear parameter s). Here, N denotes the (large) DOF of the FE space used to obtain a highly resolved solution of (25).
The main novelty here is that the proposed RB algorithm permits to compute efficiently with a reliable a posteriori error control not only the solutions of the cell problems (25) but also their derivatives with respect to the nonlinear parameter s. This is an essential ingredient to prove in Section 3.3 the uniqueness of the RB-FE-HMM macro solution and the convergence of the Newton method.
Considering an affine representation of the tensor, we first describe a suitable formulation of the cell problems before presenting the parametrized cell solution space itself. We then introduce a new a posteriori error estimator and analyze its efficiency and reliability. This is the key ingredient of the Greedy algorithm for the construction of the RB space that concludes this section.
Affine representation of the tensor. A suitable representation of the tensor where we use the transformation (11) is crucial for the RB methodology, i.e., an affine representation of the form We notice however that such direct affine representation is generally unavailable and a greedy algorithm, called the empirical interpolation method (EIM) can be used to approximate a nonaffine tensor by an affine one of the form (22) (see [14]).
Cell problems. The micro problems in the FE-HMM are based on the FE approximation of the cell functions ψ i,s K j ∈ W (K δ j ), solving the linear problem which has a unique solution using (2). For the design of the RB method, is more convenient to work in the space W (Y ) (defined in either (8) or (9)) rather than the quadrature node dependent space W (K δ j ). We thus consider the transformation (11) and using the notations the problem (23) On W (Y ) we consider the scalar product (v, w) W = Y ∇v · ∇wdy and associated norm and notice that from the ellipticity of the tensor it holds In what follows, it will be convenient to denote the micro FE space by S q (K δ j , N ) instead of S q (K δ j , T h ) to emphasize on the dimension N of the micro FE space which in RB strategy is required to be large. Analogously, the functions in S q (Y, N ) are denoted using the subscript N (e.g.,ẑ N ). The FE space S q (Y, N ) has a (shape-regular) triangulation We notice using (2) that problem (28) has a unique solution.
For the convergence of the Newton method explained in Section 3.3 we will also need to control the derivatives with respect to the parameter s of the cell functionsψ i,s K j . We assume 2 Lemma 3.1. Assume that (2) and (29) hold. Consider the solutionψ i,s N ,K j of (28). Then, where for allζ Proof. This is a standard result for FEM problems depending smoothly on a parameter (see e.g. Lemma 6.1 in [9] for details).
Parametrized cell solution space. We consider a compact subspace D of Ω × R. For any randomly chosen parameter 3 (x τ , s) ∈ D, we define the map G xτ from the physical sampling domain T δ = x τ +(−δ/2, δ/2) d centered at x τ to the reference domain Y and consider (28), (31) with a tensor a xτ ,s (y) = a ε (G xτ (y), s). Next indexed by {(T δ , s, e η ); (T δ , s) ∈ D and η = 1, · · · , d}, we define the parametrized cell solution space are the solutions of (28),(31) associated with the mapping G xτ and the Hilbert space W (Y ) is defined in either (8) or (9). On the Hilbert product space W (Y ) 2 we define the norms A posteriori error estimator. The procedure of selecting the representative cell solutions is conducted by an a posteriori error estimator which allows to control the accuracy of our output of interest (the numerically homogenized tensor) [34,19].
Assume that the RB space of dimension l, denoted by S l (Y ), is available (its construction will be detailed in Algorithm 3.4). Given the parameters (x τ , s, i), (28), (31) in S q (Y, N ) 2 and S l (Y ), respectively (i.e. with test functions (z N , ζ N ) in S q (Y, N ) 2 and S l (Y ), respectively). We then consider We derive an a posteriori estimator for bothê i,s l,T δ and ∂ sê i,s l,T δ will be analyzed in Lemma 3.3.
where the right-hand side defines a linear form on S q (Y, N ). Hence, by the Riesz theorem, there exists a uniqueē i l, We then define the residual of the a posteriori error estimator as where λ LB is an approximation of the coercivity constant λ defined in (2). We notice that the first term in (38) is the standard residual used for linear problems [33,35]. The second term arises from the nonlinearity of our problem and its control is needed to ensure the uniqueness of the nonlinear RB-FE-HMM and the convergence of the Newton method used in the implementation.
Remark 3.2. To compute the residualē i,s l,T δ in (38), we first observe that we need to solve (37), which depends on the parameter s. Thanks to the affine representation of the tensor, (37) can be decomposed into several parameter independent FE problems that can be precomputed [33,35] and hence ē i,s l,T δ W is cheap to compute. Second, for evaluating ∂ sē i,s l,T δ one can simply consider the finite difference approximation where eps is the machine precision. This can be done by solving (37) twice with parameters s and s + √ eps, respectively. In the analysis, we shall neglect the error of the above finite difference.
The next lemma gives a bound for the a posteriori error in output of interest in terms of the norms (33). It is a generalization of the result [3, Lemma 3.3] in the context of linear elliptic problems. These results are needed in our nonlinear context to control the microscopic error in the macroscopic (nonlinear) solver.
Offline algorithm. We now state step by step the offline stage of the RB algorithm. In the offline stage, we select by a greedy algorithm N triples of the form (T δn , s, η n ), where (T δn , s) belongs to a given compact D ⊂ Ω × R (since the range of the parameter s can only be obtained when the macro solution u H,RB is computed, we propose in Sect. 5 an ad hoc method to find an a priori range of s) and η n corresponds to the unit vector e ηn belonging to the canonical basis of R d . Corresponding to the N couples of (T δn , s, η n ), we computê ξ ηn,s N ,T δn , the solution of (28) with a tensor given by a xτ n ,s (y) (x τn is the barycenter of T δn ) and a right-hand side given by l ηn (·). The complete offline algorithm stated below is based on the usual procedure of the RB methodology (see [34,35]). Notice that in the case of a linear elliptic problem (i.e. a ε (x, s) independent of s), it coincides with the Greedy procedure proposed in [3].
2. Select randomly (T δ 1 , s 1 , η 1 ) ∈ Ξ RB and computeξ η 1 ,s 1 , the solution of (28) with right-hand side l η 1 (·) in S q (Y, N ), corresponding to the selected parameter (T δ 1 , s 1 , η 1 ). Set l = 1 and de- 3. For l = 2, . . . , N RB a. Compute for each (T δ , s, η) ∈ Ξ RB the residual ∆ η,s l−1,T δ defined in (38) and select the next reduced basis by choosing Remark 3.5. Notice that we choose to orthogonalize the RB in W (Y ) with respect to the scalar product (·, ·) W rather than the RB in W (Y ) 2 with respect to the scalar product (·, ·) W 2 associated to the norm · W 2 in (33) (as normally expected in the usual RB methodology) because this is more convenient in the implementation (avoiding the computation of ∂ sξ η l ,s l N ,T δ l ). Since max (T δ ,s,η)∈Ξ RB (∆ η,s l,T δ ) decays exponentially as l increases (under the assumptions of Theorem 4.3, a slight variation of the result in [20, Corollary 4.1] and [17]), we haveξ η l ,s l N ,T δ l / ∈ S l−1 and thus R l W = 0 and the above orthogonalization procedure succeeds.

Analysis of the RB-FE-HMM
In this section we first derive a priori error estimates for the RB-FE-HMM. We then show the uniqueness of the numerical approximation which is based on the convergence of the Newton method.

A priori error analysis
We introduce the following quantity to measure the error between the tensor a 0 of the homogenized problem (3) and the numerical homogenized tensor a 0 N ,K j in (17).
Assume further that (2),(4),(7) hold and that ∂ u a 0 mn ∈ W 1,∞ (Ω × R), and that the coefficients a 0 mn (x, s) are twice differentiable with respect to s, with the first and second order derivatives continuous and bounded on Ω × R, for all m, n = 1 . . . d. Then, there exist r 0 > 0 and H 0 > 0 such that, provided H ≤ H 0 and r HM M ≤ r 0 , any solution u H,RB of (13) satisfies where C is independent of H, h, N, N , ε.
Proof. We apply Lemma 7.1 (a result from [9] stated in the Appendix) withã(x K j , s) = a 0 N ,K j (s) andũ H = u H,RB .
We next have to quantify the error r HM M defined in (50) which can be decomposed as r HM M ≤ r M OD + r M IC + r RB , i.e. with the modeling, micro, and RB errors, respectively. These quantities are defined by To estimate the quantities r M IC and r M OD , we make the following smoothness and structure assumptions on the tensor: (H1) Given the degree q of the micro FE space S q (K δ j , T h ), the cell functions ψ i,s K j solution of (23) satisfy the bound |ψ i,s K j | H q+1 (K δ j ) ≤ Cε −q |K δ j |, with C independent of ε, the quadrature point x K j , the domain K δ j , and the parameter s for all i = 1 . . . d.
We next discuss the reduced basis error. Consider the space M N (Y ) as defined in (32). We want to quantify how well M N (Y ) can be approximated by the linear space S N (Y ) of dimension N . Such a quantification relies on the notion Kolmogorov N -width.
Definition 4.2. Let F be a subset of a Banach space W . We denote the distance of F to any generic N −dimensional subspace W N of W by We say that F has an exponentially small Kolmogorov N -width if there exists constants C, r > 0 independent of N such that In [20, Corollary 4.1] and [17], it is shown, for a class of symmetric linear uniformly coercive elliptic problems with continuity bound Λ and coercivity constant λ, that if the parametrized space of the RB algorithm has an exponentially small N −width (55) with constant where 0 < λ LB ≤ λ is an estimate of the coercivity bound, then the RB algorithm converges exponentially fast with respect to the dimension N of the RB space. Such assumption (55) is motivated by [31] where it is proved in the special case of one-dimensional parameter linear elliptic problems.
Notice that problem (28)-(31) with solution (ψ i,s N ,K j , ∂ sψ i,s N ,K j ) ∈ W (Y ) 2 is not coercive due to the nonlinearity of the tensor. Nevertheless, problem (28)-(31) still satisfies the following Céa inequality (see Lemma 7.3 in the Appendix) in the Hilbert space W (Y ) 2 with norm defined in (33), where we consider the solution (ψ i,s l,K j , ∂ sψ i,s l,K j ) of (28)- (31) in the space S l (Y ) (i.e. taking test functions (ẑ l ,ζ l ) ∈ S l (Y )). The above constant is given by where λ, Λ 1 , Λ 2 are the coercivity and continuity bounds (2), (29). In addition, recall the a posteriori estimate (42) of the form We obtain the following result which states that the reduced basis method converges exponentially. This is a slight adaptation of the result in [20,Corollary 4.1] and [17].
Theorem 4.3. In addition to (2) and (29), assume that the parametrized cell solution space M N in (32) has an exponentially small Komogorov N -width, with constants in (56),(58). Then, there exists constants c, κ > 0 independent of N such that for all K ∈ T H and all x K j ∈ K, whereψ i,s N ,K j andψ i,s N,K j are the solutions of the cell problem (28) in S q (Y, N ) and S N (Y ), respectively, with corresponding test functionsẑ N ∈ S q (Y, N ) Proof. Inspecting the proof of [20,Corollary 4.1] reveals that the coercivity of the problem is not needed and the Céa inequality (56) is sufficient to obtain the exponential convergence (60) of the RB algorithm using the RB space S N (Y ) constructed in the Greedy algorithm 3.4.  (13). In addition to the assumptions of Theorem 4.1, assume that u H,RB (x K j ) ∈ (u 0,low , u 0,up ), for all K ∈ T H , x K j ∈ K. Assume further (H1), (H2), and (59). Then, there exist H 0 > 0 and r 0 > 0 such at if H ≤ H 0 , h/ε ≤ r 0 , and ce −κN ≤ r 0 then for µ = 0, 1, (14), (15), (23), where r RB ≤ Λ 1 (ce −κN ) 2 with Λ 1 given in (2). We also assume δ ≤ r 0 or δ + ε/δ ≤ r 0 in the first and third cases, respectively. We use the notation H 0 (Ω) = L 2 (Ω). The constants C are independent of H, h, ε, δ, N, N .

Uniqueness of the RB-FE-HMM solution and the Newton method
Consider the derivatives with respect to s of the exact and numerical homogenized tensors in (3) and (17). We define The proof of the uniqueness of the RB-FE-HMM solution relies on the following result which is an adaptation of Lemma 4.11 in [9]. The proof of Theorem 4.5 relies on the convergence of the Newton method stated in the following lemma.
Lemma 4.6. Assume that the hypotheses of Theorem 4.5 hold. Let u H,RB be a solution of (13). Then, there exists H 0 , r 0 , ν > 0, such that provided a smallness assumption of the form (62), for all u H 0 ∈ S 0 (Ω, T H ) satisfying where C is a constant independent of H, h, k, N , N, ε. Proof. We apply Lemma 7.2 (a result from [9] stated in the Appendix) withã(x K j , s) = a 0 N ,K j (s) andũ H = u H,RB .
Proof of Theorem 4.5. The proof is an immediate consequence of Lemma 4.6, where given two numerical solutions u H,RB ,ũ H,RB of (13), we apply the Newton method with the initial guess u H 0 :=ũ H,RB . The smallness assumption (62) together with the H 1 a priori error estimate of Theorem 4.1 permits to satisfy the condition (63).
For the estimation of r HM M in (61), we consider the decomposition r HM M ≤ r M OD + r M IC + r RB where r M OD , r M IC , r RB are defined similarly to (52),(53),(54), respectively, with the exception that all the tensors are differentiated with respect to the s parameter. Hence hypothesis (H1') and (H2'), similar to (H1) and (H2) but for ∂ s ψ i,s K j are needed. Following [8], the quantities r M OD , r M IC satisfy analogous estimates to those of r M OD , r M IC . It remains to estimate This is done in the following lemma.
Lemma 4.7. Assume that the hypotheses of Theorem 4.4 hold with a periodic coupling with δ = ε. Assume further (29), and (59). Then, there exist constants c, κ > 0 such that Proof. Differentiating (46) with respect to s and using Theorem 4.3 conclude the proof.
Remark 4.8. Notice that a similar a posteriori estimator as used here in the offline stage to control r RB and r RB could be used to define a RB-FE-HMM for linear parabolic multiscale problems with a time-dependent tensor of the form ∂u ε ∂t (x, t) = ∇ · (a ε (x, t)∇u ε (x, t)) + f (x, t) as analysed in [10]. In this case, the micro problems would be parametrized by the location of the cell problem in the domain Ω and the time variable of the tensor a ε .

Implementation issues
We give in this section details on an efficient implementation of the offline and online procedures of the proposed RB-FE-HMM.

Offline procedure
Output of the offline procedure. The output of the offline procedure is the RB space (10). Rather than storing the reduced basis functions, using the affine representation (22) described above, it is sufficient to compute the following matrices and vectors We emphasize that offline stage is only operated once and the output of offline stage can be repeatedly used for the online stage computations independently of the chosen online solver.
Preprocess of the offline stage. In Algorithm 3.4, we assume that the value of the solution u H (x τ ) at the quadrature point x τ lies in a bounded real interval for each τ in training set. But the range of u H can only be known once the macro solution is computed due to the nonlinearity of the problem. One possibility to address this issue is to map R into a bounded interval. However, our numerical tests indicate that this approach fails in general unless an unreasonably large training set is used. Therefore, we propose to use an empirical algorithm to get an a priori guess of the solution range, motivated by the following classical result [29, Chapter 1.6].
We then apply the following procedure: we first solve (3) replacing the homogenized tensor a 0 (x τ , s) alternatively with Y a xτ ,s (y) −1 dy −1 and Y a xτ ,s (y)dy. We then set a maximum range (u 0,low , u 0,up ) by taking the minimum and maximum values of the two solutions obtained from the first step. In the end, we consider the range of parameter s to be (u 0,low −α, u 0,up +α) enlarged by ±α (a safety factor of about 10%).

Online procedure
We describe here how the online stage of the RB-FE-HMM can be efficiently implemented.
Fast solution of micro problems. Owing to the affine form (22) of the tensor a ε , the problem (15) amounts to solving an N × N linear system (recall that N , the number of RB for all m = 1, . . . N. Next, again thanks to the affine representation of the tensor (here we are assuming the representation (22) for simplicity), (66) can be written as or equivalently where the N × N matrices A p , p = 1, . . . , P and the vectors F i p ∈ R N , p = 1, . . . , P, i = 1, . . . , d are defined by (64).
We emphasize that the matrices A p and the vectors F i p are assembled and stored in the offline stage, thus (68) amounts just in building the linear combination by evaluating Θ p (·, ·) at the desired parameter (x K j , u H (x K j )) for the tensor a xτ ,s (y) and solving the N ×N system (68) for each micro function at the quadrature points needed to assemble (14).
Newton method implementation. We consider a sequence of {u H k } in S 0 (Ω, T H ) and express each function in the FE basis of S 0 (Ω, The Newton method (20) can be written out in terms of matrices as where B(u H k ), B (u H k ) are the stiffness matrices associated to the bilinear forms B H (z H ; ·, ·), B H (z H ; ·, ·) defined in (14) and (19), respectively. Here, F is a vector associated to the source term (13), which also contains the boundary data if for instance general Dirichlet or Neumann boundary conditions are considered.
Stiffness matrices. Following the implementation in [6] we consider for each element K ∈ T H the FE basis functions {φ H K,i } n K i=1 associated with this element and the local contribution where ϕ are the solutions of (15) constrained by φ H K,p , φ H K,q , linearized at x K j , respectively.
Similarly to the FE-HMM implementation [9], we see by differentiating (70) that the stiffness matrix B (U ) in (69) associated to the non-symmetric form B H (z H ; ·, ·) defined in (19) is given by the sum of J products of n K × n K matrices where the column vector U K,k of size n K gives the components of z H in the basis of the macro element K ∈ T H . Here, the derivative with respect to s of the n K × n K matrix B K,j (s) can be simply approximated by the finite difference ∂ ∂s (B K,j (s)) ≈ , where eps is the machine precision. Therefore, the cost of computing the stiffness matrices for both B(u H k ) and B (u H k ) is about twice the cost of computing the stiffness matrix B(u H k ) alone.

Numerical experiments
This section is dedicated to the numerical illustration of the theoretical results in Sect. 4 and the performances of the proposed RB-FE-HMM compared to the FE-HMM. All computations are performed in Matlab on the same machine (with a non-parallel implementation). We first consider a simple illustrative example and then the stationary Richards problem which has a nonaffine tensor.

Numerical evaluation of the errors.
Let u H be the numerical solution and u ref be a reference solution (for the effective problem (1)) computed on a fine triangulation T h of Ω. The errors u ref − u H in the H 1 and L 2 norms are estimated by where {z K j , ρ K j } denotes the quadrature formula on the fine triangulation T h .
Macro FEM and QF used in the examples. In the following examples, when using P1 triangular (tetrahedral) elements for the macro problems, we choose the barycenter of the element as single quadrature point and the weightω = |K|. When we use P2 triangular elements for the macro problems, we choose the Gauss three points quadrature formula with barycentric coordinates (1/6, 1/6, 2/3) and weightsω i = |K|/3, i = 1, 2, 3.

A simple illustrative example
We consider the model problem (1) in Ω = [0, 1] 2 with f = 50e (x 1 −0.2) 2 +(x 2 −0.3) 2 , and the following mixed boundary conditions, Consider a diagonal multiscale tensor with the following affine expression 6 a ε (x, s) 11 = (x 2 1 + 0.2) + (x 2 sin(sπ) + 2)(sin(2π Using the homogenization theory [29], the corresponding homogenized tensor is also diagonal with entries given by the harmonic averages Offline stage. In the offline stage, we set the parameter space to be D = Ω i × U , where Ω i is a closed subset of Ω such thatT δ = x τ + [−δ/2, δ/2] d ⊂Ω for all τ ∈ Ω i , and U is a closed bounded interval of R (an estimation of the range for u 0 ). In order to obtain U , motivated by Lemma 5.1, we first solve (1) on a coarse 8 × 8 macro mesh by replacing the homogenized tensor respectively with the arithmetic and harmonic averages of the multiscale tensor. The ranges of the corresponding solutions are shown in Table 1 and we choose U = [0.9, 3.66] adding a safety correction. For the RB offline stage, we propose in Sect. 3.4 a new a posteriori error estimator (38) in order to guarantee the convergence of the Newton method. We will also check the computational overhead of this new estimator compared to the (standard) a posteriori estimator used for linear problems [3]. The offline parameters are collected in Table 2 and the comparison of the two estimators are shown in Table 3. We observe in this test that the offline output using ∆ l,T δ has only one additional basis function.
Online stage: convergence rates for the P1 and P2 RB-FE-HMM. Using the computed offline output (obtained by via ∆ l,T δ as the offline estimator), we consider a P1 FEM and a P2 FEM for the online stage. The reference solution u ref ≈ u 0 is obtained using the P2 FEM with a 1024 × 1024 uniform mesh.   By the a priori estimates of Theorems 4.1-4.4 and the mesh size used for the offline computation, we have the bound ∆ 2 l,T δ = O(tol RB ) ∼ 10 −10 for r RB and r M IC ∼ 10 −7 . As we choose sampling domains with size δ = ε with periodic boundary condition we have r M OD = 0 and we expect r HM M ∼ 10 −7 . We observe in Fig. 1 the expected convergence rates with respect to the macro mesh.
RB-FE-HMM v.s FE-HMM. In this test, we compare the efficiency and accuracy between the P1 RB-FE-HMM and the P1 FE-HMM. For the P1 FE-HMM, we set N M IC = N M AC for each refinement step (L 2 refinement strategy), where N M AC and N M IC are the numbers of macro and micro DOF in one space direction, respectively. We can see in Table 4 that the H 1 and L 2 errors for the two methods decay with the same rates which is consistent with the a priori estimates. However, the RB-FE-HMM has a considerably reduced computational cost for fine meshes (up to two orders of magnitude in this example). Next we present = 1300s and we denote the total offline and online time by t total RB . We see that the FE-HMM is still significantly more expensive except for coarse macroscopic meshes. This indicates that even for one computation, the RB-FE-HMM can provide an important computational speed-up.

Stationary Richards problem
We consider the stationary Richards equation for describing the fluid pressure in an unsaturated porous media with a nonlinear permeability tensor K ε (s) similar to the one in [21, Sect. 5.1] written as .
Notice that this problem can be cast in the form (1) by using the change of variable v ε (x) = u ε (x)−x 2 . We set f = 1 and consider Dirichlet conditions on the top boundary of the domain and Neumann conditions on the rest of boundaries, that is Offline stage. The parameters are given in Table 6. Similarly as explained for the previous example, we determine an a priori range for the homogenized solution U = [0.9, 3.93]. As the permeability tensor K ε does not have an affine representation (22), we need to apply the EIM, which introduces another error term r EIM in r HM M , as discussed in [3] (this term can be controlled by the prescribed tolerance tol EIM [14]).
Online stage. We plot in Fig. 2 the online solution u H,RB on a uniform 65 × 65 macro mesh for the P1 RB-FE-HMM (left picture) and the FE-HMM (right picture). We observe that the two solutions are very similar as expected.
We notice that the range of u H,RB is [−2.9, −1] which safely lies in our a priori range [−3.1, −0.8]. Therefore, the offline output can be successfully used for the online computation. In Table 7, we present the errors of the Newton method iterations and the corresponding CPU times. Due to the quadratic convergence rate of the Newton method, only four iterations are needed in all considered cases to reach the machine precision.

Appendix
The proof of Theorem 4.1 relies on the following lemma taken from [9] and based on the analysis for standard FEM with numerical quadrature from [11]. It is a reformulation of the statement of Theorem 3.1 in [9], its proof is thus omitted. ω K,jã (x K j ,ũ H (x K j ))∇ũ H (x K j )·∇w H (x K j ) = Ω f (x)w H (x)dx, ∀w H ∈ S 0 (Ω, T H ), we have the H 1 and L 2 error estimates where C is independent of H, Q H and the tensorã.
The proof of Theorem 4.5 relies on the following result which is a reformulation of Lemma 4.11 in [9]. Lemma 7.2. Assume that the hypotheses of Lemma 7.1 hold. Assume further thatã(x, s) is twice continuously differentiable with respect to s with derivatives continuous and bounded on Ω × R. Then, there exists H 0 , R 0 , ν > 0, such that for Finally, for the proof of Theorem 4.3, we need to prove the Céa inequality (56) for the problem (28)-(31). Lemma 7.3. Assume (2) and (29). Then (56) holds with constant C 0 in (57).
Proof. We denoteê =ψ i,s N ,K j −ψ i,s l,K j and ∂ sê = ∂ sψ i,s N ,K j − ∂ sψ i,s l,K j . Considering the problem (28)- (31) with test functions in S q (Y, N ) 2 and S l (Y ), respectively, and subtracting, we deduce for all (z l , ζ l ) ∈ S l , b(∂ sê ,ζ l ) = − Y ∂ s a x K j ,s (y)∇ê(y) · ∇ζ l dy, where the symmetric bilinear form b(·, ·) is defined in (24). Using We deduce from the Young inequality and (2), Using again the Young inequality yields Finally, the application of the Céa lemma to (28) yields ê W ≤ Λ 1 λ inf z l ∈S l (Y ) ψ i,s N ,K j −z l W which permits to conclude the proof.