Error estimate of the Non Intrusive Reduced Basis method with finite volume schemes

The context of this paper is the simulation of parameter-dependent partial differential equations (PDEs). When the aim is to solve such PDEs for a large number of parameter values, Reduced Basis Methods (RBM) are often used to reduce computational costs of a classical high fidelity code based on Finite Element Method (FEM), Finite Volume (FVM) or Spectral methods. The efficient implementation of most of these RBM requires to modify this high fidelity code, which cannot be done, for example in an industrial context if the high fidelity code is only accessible as a"black-box"solver. The Non Intrusive Reduced Basis method (NIRB) has been introduced in the context of finite elements as a good alternative to reduce the implementation costs of these parameter-dependent problems. The method is efficient in other contexts than the FEM one, like with finite volume schemes, which are more often used in an industrial environment. In this case, some adaptations need to be done as the degrees of freedom in FV methods have different meenings. At this time, error estimates have only been studied with FEM solvers. In this paper, we present a generalisation of the NIRB method to Finite Volume schemes and we show that estimates established for FEM solvers also hold in the FVM setting. We first prove our results for the hybrid-Mimetic Finite Difference method (hMFD), which is part the Hybrid Mixed Mimetic methods (HMM) family. Then, we explain how these results apply more generally to other FV schemes. Some of them are specified, such as the Two Point Flux Approximation (TPFA).


The non-intrusive reduced basis method
Let Ω be an open bounded domain in R with ≤ 3. The NIRB method in the context of a high fidelity solver of finite element or finite volume types involves two partitioned meshes, one fine mesh ℳ ℎ and one coarse mesh ℳ , where ℎ and are the respective sizes of the meshes and ℎ ≪ . The size ℎ (res. ) is defined as ℎ = max ∈ℳ ℎ ℎ (resp. = max ∈ℳ ) where the diameter ℎ (or ) of any element in a mesh is equal to sup , ∈ | − |.
As it is classical in other reduced basis methods, the NIRB method is based on the assumption (assumed or actually checked) that the manifold of all solutions = { ( ), ∈ } has a small Kolmogorov -width ( ) [24]. This leads to the fact that very few well chosen solutions are sufficient to approximate well any element in . These well chosen elements are called the snapshots. In this frame, the method is based on two steps (see Fig. 1): one offline step and one online. The "offline" part is costly in time because the snapshots must be generated with a high fidelity code on the fine mesh ℳ ℎ . The "online" step is performed on the coarse mesh ℳ , and thus much less expensive than a high fidelity computation. This algorithm remains effective as the offline part is performed only once and in advance and also independently from the online stage. The online stage can then be done as many times as desired.
-In the offline part, several snapshots are computed on the fine mesh for different well chosen parameters in the parameter set with the (fine and costly) solver. The best way to determine the required parameters is through a greedy procedure [1,6,30] if available or through an Singular Value Decomposition (SVD) approach. -The online part consists in computing a coarse solution with the same solver for some (new) parameter ∈ and then 2 -project this (coarse) solution on the (fine) reduced basis. This results in an improved approximation, in the sense that we may retrieve almost fine error estimates with a much lower computational cost. This is not as in classical extrapolation schemes since the solution on the fine mesh is not employed to obtain the approximation.

Motivation and earlier works
Several papers have underlined the efficiency of the NIRB method in the finite element context, illustrated both with numerical results presenting error plots and the online part computational time [8][9][10]25]. However, to the best of our knowledge, works with Finite Volume (FV) schemes have not yet been studied with a nonintrusive approach [7,21,23,28,29], and they are often preferred to finite element methods in an industrial context. Thanks to recent works on super-convergence [16], and with some technical subtleties, we are now able to generalize the two-grids method which is non-intrusive to FV methods and propose the numerical analysis of this method.

Main results
In the context of 1 -FEM solvers, the works [8,25] retrieve an estimate error of the order of (ℎ + 2 ) in the energy norm using the Aubin-Nitsche's lemma [3] for the coarse grids solution (for a reduced basis dimension large enough). With FV schemes, no equivalent of the Aubin-Nitsche's lemma is available, instead, we consider the class of Hybrid Mimetic Mixed (HMM) schemes for elliptic equations and use a super-convergence property proven in [13,14,16]. We will first focus on hMFD scheme, which is part of the family of Hybrid Mimetic Mixed methods (HMM) [12,13,[17][18][19]. It is a finite volume method despite its name. Indeed hMFD scheme relies on both a flux balance equation and on a local conservativity of numerical fluxes. It uses interface values and fluxes as unknowns. The particularity of hMFD is that the cell unknowns must be located at the center of mass of the cells [4,5]. HMM also includes mixed finite volume schemes (MFV) [15] and hybrid finite volume schemes (HFV), a hybrid version of the SUSHI scheme [20]. All HMM are built on a general mesh, namely a polytopal mesh, which is a star-shaped mesh regarding the unknowns of the cells.
Let us consider the following linear second-order parameter dependent problem as our model problem: where ∈ 2 (Ω), is a parameter in a set , and for any ∈ , (.; ) : Ω → R × is measurable, bounded, uniformly elliptic, and (x; ) is symmetric for a.e. x ∈ Ω. Under general hypotheses, it is well-known that this problem has a unique solution.
The usual weak formulation for problems (1.1a), (1.1b) reads: Find ∈ 1 0 (Ω) such that, ∀ ∈ 1 0 (Ω), ( , ; ) = ( , ), The main result of this paper is the following estimate: Theorem 1.1 (NIRB error estimate for hMFD solvers). Let ℎ ( ) be the reduced solution projected on the fine mesh and generated with the hMFD solver with the unknowns defined on x = x (the cell centers of mass), and ( ) be the exact solution of (1.2) under an 2 regularity assumption (2.5) (which will be stated later), then the following estimate holds where 1 and 2 are constants independent of ℎ and , 2 depends on , the number of functions in the basis, and ‖·‖ is the discrete norm introduced in Section 2, and depends of the Kolmogorov -width. If is such as 2 ∼ ℎ, and ( ) small enough, it results in an error estimate in (ℎ).
Note that if is chosen such as 2 ∼ ℎ and ( ) small enough, it results in an error estimate in (ℎ).

Outline of the paper
The rest of this paper is organized as follows. In Section 2, we describe the mathematical context. In Section 3, we recall the two-grids method. Section 4 is devoted to the proof of Theorem 1.1 with the hybrid-Mimetic Finite Difference scheme (hMFD). Section 5 generalizes Theorem 1.1 to other schemes, such as the Two Point Flux Approximation (TPFA). In the last section, the implementation is discussed and we illustrate the estimate with several numerical results on the NIRB method.

Mathematical background
In this section, we recall the definition of the Hybrid Mixed Mimetic methods (HMM) family, some notations and some of its properties [12,18,19] that will be necessary for the analysis of NIRB method in this finite volume context.

The Hybrid Mixed Mimetic methods (HMM) family
Describing the HMM family requires to introduce the Gradient Discretisation (GD) method [19], which is a general framework for the definition and the convergence analysis of many numerical methods (finite element, finite volume, mimetic finite difference methods, etc.).
The GD schemes involve a discete space, a reconstruction operator and a gradient operator, which taken together are called a Gradient Discretisation. Selecting the gradient discretisation mostly depends on the boundary conditions (BCs). We now introduce the definition of GD for Dirichlet BCs as in [19] and the GD scheme associated to our model problem. In what follows, we will refer to Π or Π ℎ depending on the mesh considered and for the gradient reconstruction too (respectively ∇ or ∇ ℎ ).
We will use two general polytopal meshes ( [19], Def. 7.2) which are admissible meshes for the hMFD scheme. (1) ℳ is a finite family of non-empty connected polytopal open disjoint subsets (the cells) such that Ω = ∪ ∈ℳ . For any ∈ ℳ, | | > 0 is the measure of and ℎ denotes the diameter of . (2) ℱ = ℱ int ∪ℱ ext is a finite family of disjoint subsets of Ω (the edges of the mesh in 2D), such that any ∈ ℱ int is contained in Ω and any ∈ ℱ ext is contained in Ω. Each ∈ ℱ is assumed to be a nonempty open subset of a hyperplane of R , with a positive ( − 1)-dimensional measure | |. Furthermore, for all ∈ ℳ, there exists a subset ℱ of ℱ such that = ∪ ∈ℱ . We assume that for all ∈ ℱ, ℳ = { ∈ ℳ : ∈ ℱ } has exactly one element and ⊂ Ω or ℳ has two elements and ⊂ Ω. The center of mass is x , and, for ∈ ℳ and ∈ ℱ , n , is the (constant) unit vector normal to outward to .
is a family of points of Ω indexed by ℳ and ℱ, denoted by = ((x ) ∈ℳ , (x ) ∈ℱ ), such that for all ∈ ℳ, x ∈ and for all ∈ ℱ, x ∈ . We then denote by , the signed orthogonal distance between x and ∈ ℱ , that is: for all x ∈ . We then assume that each cell ∈ ℳ is strictly star-shaped with respect to x , that is , > 0 for all ∈ ℱ . This implies that for all x ∈ , the line segment [x , x] is included in . We denote x the center of mass of and by x the one of . For all ∈ ℳ and ∈ ℱ , we denote by , the cone with vertex x and basis , that is (4) is a set of points (the vertices of the mesh). For ∈ ℳ, the set of vertices of , i.e. the vertices contained in , is denoted . Similarly, the set of vertices of ∈ is .
The Figure 2 illustrates a cell of a 2D polytopal mesh. The regularity factor for the mesh is In what follows, we will consider two polytopal meshes. The fine mesh will be denoted ℎ = (ℳ ℎ , ℱ ℎ , ℎ , ℎ ) and = (ℳ , ℱ , , ) will be referred to as the coarse mesh. All HMM schemes require to choose one point inside each mesh cell x , and in the case the center of mass x is chosen, then the scheme corresponds to hMFD and superconvergence is well known [12,13,16]. Until Section 5, we will consider x = x . Definition 2.4 (Hybrid Mimetic Mixed gradient discretisation (HMM-GD)). For hMFD scheme, we use the following GD ( [19], Def. 13.1.1): Π : ,0 → 2 (Ω) is the following piecewise constant reconstruction on the mesh: As explained in the introduction of this chapter, hMFD, HFV and MFV schemes are three different presentations of the same method. With the notations above, any HMM method for the weak form (1.2) can be written ( [17], Eq. (2.25)): Find where is our variable parameter, ( ) is the 2 projection of ( ) on and B = ((B ) , ′ ) , ′ ∈ℱ is a symmetric positive definite matrix, resulting from the definition of ∇ .

The hybrid Mimetic Finite Difference (hMFD) method
We now introduce the super-convergence property on hMFD which will be used in the proof of Theorem 1.1, but first we need the following 2 regularity assumption (which holds if is Lipschitz continuous and Ω is convex): Let ∈ 2 (Ω), the solution ( ) to (1.2) belongs to 2 (Ω), and with depending only on Ω and . We define ℳ ℎ : 2 (Ω) → 2 (Ω) as the orthogonal projection on the piecewise constant functions on ℳ ℎ that is Theorem 2.5 (Super-convergence for hMFD schemes ( [16], Thm. 4.7)). Let ≤ 3, ∈ 1 (Ω), and ( ) be the solution of (1.2) under assumption (2.5). Let ℎ be a polytopal mesh, and be an HMM gradient discretisation on ℎ with the unknowns defined on x , and let ℎ ( ) be the solution of the corresponding GD. Recall that x is the center of mass of and we are in the case where x = x . Then, considering ( ) as the piecewise constant function on ℳ ℎ equal to (x ; ) on ∈ ℳ, there exists > 0 not depending on ℎ such that To recover (2.6) in the case x = x , we used the Lemma 7.5 of [16] on the approximation of 2 functions by affine functions to obtain Remark 2.6. We consider here ‖·‖ as the discrete semi norm of 1 so as not to make notations too cumbersome. The usual discrete semi-norm for 1 is defined by Under some conditions on the regularity of the mesh, this norm and ‖∇ ·‖ 2 (Ω) are equivalent ( [19], Lem. 13.11).
In the next section, we recall the offline and the online parts of the two-grids algorithm.

Main steps
This section recalls the main steps of the two-grids method algorithm [8,25]. Let ℎ ( ) refer to the hMFD solution on a fine polytopal mesh ℎ , with cells ℳ ℎ and respectively ( ) the one on a coarse mesh , with the cells ℳ . We briefly recall the NIRB method. Points 1 and 2 are performed in the offline part, and the others are done online.

Rectification post-process
We introduce ℎ ( ) = ∫︀ Ω Π ℎ ℎ ( ) · Π ℎ Φ ℎ dx. The rectification process, explained in [8,10,25], can be employed in addition of the NIRB classical algorithm. This implies that if the true solution is in the reduced space, then the NIRB method will give this true solution. Let be the number of parameters in . Let A be the matrix such that , = ( ), ∀ ∈ , ∀ = 1, . . . , , and B be the matrix such that , = ℎ ( ), ∀ ∈ , ∀ = 1, . . . , . The aim is to minimize, on R , ‖AR − B ‖ 2 . The solution of this problem with a regularization parameter is the rectification matrix: where is the regularization parameter.
The approximation used with this post-process is In the next section, we detail how to obtain the classical finite elements estimate in (ℎ) on the NIRB algorithm, when the snapshots are computed with the hMFD GD using a polytopal mesh.

NIRB error estimate
In this section, we consider x = x which is the case with the hMFD scheme. Some other cases will be detailed in Section 5.
We now continue with the proof of Theorem 1.1.
Proof. In this proof, we will denote for ≤ with not depending on ℎ or . We use the triangle inequality on -The first term 1 can be estimated using a classical result for finite volume schemes (consequence of [19], Prop. 13.16) such that: -The best achievable error in the uniform sense of a fine solution projected into ℎ relies on the notion of Kolmogorov -width ( [26], Thm. 20.1). If is a compact set in a Banach space , the Kolmogorov -width of is Here we suppose the set of all the reconstructions of the solutions = {Π ℎ ℎ ( ), ∈ } has a low complexity which means for an accuracy = ( ) related to the Kolmogorov -width of the manifold , there exists a set of parameters { 1 , . . . , } ∈ , such that [6,8,11,25] -Consider the term 3 now. We will need the following proposition where the property of super-convergence for the hMFD scheme (2.6) is used.
Proof. Since ℳ is a partition of Ω, To begin with, let Π 0 : (Ω) → ∞ (Ω) be the piecewise constant projection operator on ℳ such that: We use the triangle inequality on the left part of the inequality (4.5) and therefore, (4.8) -We first consider the term 3,1 . But first, this requires the use of a further operator which we now introduce. Each cell ∈ ℳ is star-shaped with respect to a ball centered in x of radius = min ∈ℱ , ( [19], Lem. B.1). We then use an averaged Taylor polynomial as in [3] but simplified. Let us consider the following polynomial of ( ) averaged over : This polynomial is of degree less or equal to 1 in x. Let us introduce Π 1 : 1 (Ω) ∩ (Ω) → R, the piecewise affine projection operator on ℳ such that: With the triangle inequality, we obtain (4.11) • Using the Cauchy-Schwarz inequality, (4.14) Thus, with the inequalities (4.14) and (4.13), we get taking the square and integrating over , we obtain 16) and summing over yields The inequality (4.17), combined with (4.12), entails that • The term 3,1,2 can be estimated using a continuous reconstruction of Φ ℎ , denoted by Φ .

Results on other FV schemes
In this section, we consider the case where x is not the center of mass, as it is the case for some FV schemes. Therefore the left hand side of the inequality (4.22) cannot be estimated using equation (4.20). The unknowns x are not necessarily the centers of mass of the cells neither with HMM methods nor with the TPFA scheme [2,14]. Under the following superadmissibility condition the TPFA scheme is a member of the the HMM family schemes ( [19], Sect. 13.3, [17], Sect. 5.3) with the choice ℒ = . This leads to take x as the circumcenters of the cells with 2D triangular meshes. Theorem 1.1 holds in 2D on uniform rectangles with TPFA since the superadmissibility condition is satisfied in this case where x is the centre of mass of the cells. The TPFA scheme is rather simple to implement, and therefore we will present in the last section numerical results with a TPFA solver. We will use the definition of a local grouping of the cells as in [16] (Def. 5.1). We will extend the Theorem 1.1 in the case where such groupings of cells exist.
Definition 5.1. (Local grouping of the cells). Let be a polytopal mesh of Ω. A local grouping of the cells of is a partition G of ℳ , such that for each ∈ G, letting := ∪ ∈ , there exists a ball ⊂ such that is star-shaped with respect to . This implies that for all x ∈ and all y ∈ , the line segment [x, y] is included in . We then define the regularity factor of G and, with = x − x , and Note that we are interested in situations where |e | = ⃒ is much smaller than |e | ∀ ∈ . The aim of this section is to estimate the left hand side of the inequality (4.22) in ( 2 ) using a local grouping of the cells. The rest of the proof remains unchanged.
We will need the following theorem of super-convergence for HMM schemes with local grouping ( [16], Thm. 5.4).
Theorem 5.2 (Super-convergence for HMM schemes with local grouping ( [16], Thm. 5.4)). Let ∈ 1 (Ω), and ( ) be the solution of (1.2) under assumption (2.5). Let ℎ be a polytopal mesh, and be an HMM gradient discretisation on ℎ and G be a local grouping, and let ℎ ( ) be the solution of the corresponding GD. Then, considering ( ) as the piecewise constant function on ℳ ℎ equal to (x ; ) on ∈ ℳ, there exists not depending on or ℎ such that Theorem 5.3 (NIRB error estimate with local grouping). Let ℎ ( ) be the reduced solution projected on the fine mesh and generated with the hMFD solver with the unknowns defined on x such that G is in ( 2 ) on the coarse mesh, and ( ) be the exact solution of (1.2) under assumption (2.5), then the following estimate holds where 1 and 2 are constants independent of ℎ and , 2 depends on , the number of functions in the basis, and ‖·‖ is the discrete norm introduced in Section 2, and depends of the Kolmogorov -width. If is such as 2 ∼ ℎ, and ( ) small enough, it results in an error estimate in (ℎ).
Proof. In this proof, we will still denote for ≤ with not depending on ℎ or . The reconstruction Φ of Φ ℎ must belong to 1,∞ . As in the previous section, with the equation (4.20), As in the previous section (4.23), Thus, the inequality (5.7) yields With the triangle inequality, the first term in (5.9) becomes Using the decomposition of the mesh in patches and with the definition of , the first term of (5.10) gives (5.11) . This concludes the proof since the rest is similar to the one of Theorem 1.1. Note that for the estimate of 3,2 (4.30), the equation (5.5) from the Theorem of super-convergence with local grouping is used instead of (2.6).

Some details on the implementation and numerical results
We consider two simple cases in 2D for the numerical results with the TPFA scheme. Both results are computed on the unit square. We use an harmonic averaging of the diffusion coefficient ( [17], Sect. 5.3). Our variable parameter is ∈ R 4 = ( 1 , 2 , 3 , 4 ). For both cases, the size of the meshes is defined as the maximum length of the edges. The diffusion coefficient we consider here is ( ) = (2 1 + 2 sin( + ) cos( )) and = ( 3 (1 − ) + 4 (1 − )). We choose random coefficients in [0, 1] for the snapshots with = 5 and our solution is defined with 1 = 0.99, 2 = 0.8, 3 = 0.2, 4 = 0.78. For the exact solution, we consider the TPFA solution on a finer mesh (Figs. 3 and 4). For the computation of the norm, we use the discrete semi-norm as in the remark of the Section 2 (2.7). NIRB results (with and without the rectification 3.2) are compared to the classical FV errors (Figs. 5 and 6). We measure the following relative error In practice, one approach, based on the computation times, consists in choosing a precise time 1 and in finding the associated coarse solution computed within this time. Then, the fine grid is chosen such that 2 = ℎ. In our tests, we choose several fine mesh sizes to analyze the rate of the error, and the coarse mesh size is equal to 0.25. An other approach can be to select a fine mesh size such that the method works for several coarse mesh sizes.

Uniform grid
The first case presents results on a rectangular uniform grid where x is the center of mass of the cell.

Triangular mesh
The second case is defined on a triangular mesh where x are the circumcenter of the cells, such that G is in ( 2 ).

Discussion on the implementation
We implemented the TPFA scheme on Scilab and retrieved several solutions for the NIRB algorithm on Python to highlight the black box side of the solver. The scilab files consist in three text files with solution values, the cell center coordinates, and one file with information on the edges (distance , the area between the cell center and the edge, and the labels).
-Implementation of the TPFA method. , on ℱ ext . To assemble the matrices of the TPFA scheme, we iterate on each edge, and we add the harmonic average on each cell, and for we add the term | , | × ( ).
-Time execution (min, s). Remark 6.1. In dimension 2, we expect a speedup of 1/ℎ. Indeed, the degrees of freedom ℎ (for the fine mesh) are of order (1/ℎ) 2 (resp. = (1/ ) 2 for the coarse mesh), and the costs of an optimal solver are in ( ℎ ) (or ( ) for the coarse mesh). Thus the speedup with ℎ = 2 is equal to 1/ℎ and differs from other classical reduced-basis methods. In our case, this is difficult to observe since our model problem is very simple with few degrees of freedom, and the computational costs take into account other subroutines such as mesh readers which are not proportional. Remark 6.2. Note that for the discontinuous diffusion coefficient , with the TPFA scheme, we recovered numerically the same estimate as in the Lipschitz continuous case, when we use the harmonic mean even if the proof no longer works.