A posteriori error estimation for numerical model reduction in computational homogenization of porous media

Numerical model reduction is adopted for solving the microscale problem that arizes from computational homogenization of a model problem of porous media with displacement and pressure as unknown fields. A reduced basis is obtained for the pressure field using (i) spectral decomposition (SD) and (ii) proper orthogonal decomposition (POD). This strategy has been used in previous work—the main contribution of this article is the extension with an a posteriori estimator for assessing the error in (i) energy norm and in (ii) a given quantity of interest. The error estimator builds on previous work by the authors; the novelty presented in this article is the generalization of the estimator to a coupled problem, and, more importantly, to accommodate the estimator for a POD basis rather than the SD basis. Guaranteed, fully computable and low‐cost bounds are derived and the performance of the error estimates is demonstrated via numerical results.

solution space of a discrete RVE problem. In particular, we highlight methods based on superposition of "modes" that are characteristic to the solution field. Various attempts have been made, in the context of small strain (visco)plasticity, to approximate the inelastic strains by "inelastic modes." One example is the "eigendeformation-based-reduced-order homogenization" technique introduced by Fish and coworkers, 3,4 which relies on the concept of transformation field analysis, originally proposed by Dvorak and Benveniste. 5 Michel and Suquet 6,7 proposed a similar approach coined nonuniform transformation field analysis (NTFA). Fritzen et al [8][9][10][11] exploited NTFA combined with proper orthogonal decomposition (POD) for viscoelasticity and a class of standard dissipative materials. Jänicke et al 12 applied this approach to poroelasticity, whereby the pore pressure plays a role similar to inelastic strains in the NTFA framework. As a consequence of the reduced model, the macroscale problem reduces to a single-phase continuum, where the "mode coefficients" can be interpreted as internal variables. This is the reduction model that is considered in this article. Other reduction techniques have been studied in the context of multiscale finite element methods by, e.g. Nguyen 13 for parametrized PDEs and by Efendiev et al 14,15 for flows in heterogeneous media. The concept of "hyperreduction" has also been investigated for multiscale methods by, e.g. Hernández et al 16 and Memarnahavandi et al,17 where, in addition to reducing of the number of degrees of freedom, also the cost for evaluating the residual is reduced.
Naturally the use of a reduced basis results in an additional source of error, and therefore it is of interest to control and quantify this error. Several different error estimators have been developed, for different model reduction techniques, in the context of multiscale modeling. Methods for estimating the error from POD type reduction techniques have been presented by Abdulle et al 18,19 for heterogeneous multiscale methods and by Boyaval 20 for numerical homogenization. Ohlberger and Schindler developed a method for estimating the error for localized reduced basis multiscale methods. Error estimation based on the constitutive relation error has been proposed by Kerfriden et al 21 and Chamoin and Legoll. 22 For the estimator proposed in this article, we consider the fully resolved finite element solution to be the exact one, similar to what was presented in Aggestam et al 23 and Ekre et al. 24 The consequence of this is that the estimator only quantifies the error stemming from the use of a reduced base and "ignores" e.g. space and time discretization errors. Naturally it would be of interest to develop an estimator which quantifies these errors as well. However, there are cases where the point of departure is already a fine mesh; it might, for example, be obtained from detailed voxel data or be required to capture complex microstructural features. In addition, for linear problems, many of the necessary quantities can be precomputed in an "offline" stage where one can afford the fine mesh without having to tackle the problem of integration, as one would have to do for, e.g. a nonlinear problem.
In this article, we consider a continuum mechanics model of porous media as the model problem, following the work in Jänicke et al. 12 We will consider two different methods of obtaining a reduced basis for the microscale problem; spectral decomposition (SD), which was used for the linear transient heat flow problem in, e.g., Aggestam et al, 23 Ekre et al, 24 and POD, which was used for the same porous media problem in Reference 12. The main new contribution of this article, as compared with Jänicke et al, 12 is the extension with an a posteriori error estimator for assessing the accuracy of the reduced solution. We utilize the linearity of the problem and derive an auxiliary symmetric form, cf. e.g. Pares et al. [25][26][27] From the auxiliary form we obtain error bounds based on the discrete residual, cf. e.g. Jakobsson et al. 28 We aim for guaranteed, explicit, fully computable, bounds on the error introduced by the reduced basis in terms of (i) energy norm and (ii) a user-defined quantity of interest using goal oriented error estimation, cf. Oden and Prudhomme. 29 We aim to base the estimate on the "active" modes, that is, the modes used for the reduced solution, and will thus not consider hierarchical approaches (eg, by computing additional modes used solely for the purpose of estimating the error). This strategy for estimating the error follows what was presented in References 23,24. The novelty in this article is the generalization of the estimator. They key differences can be summarized as follows; the estimator is derived for a POD basis (rather than for a SD basis), the estimator is generalized for the coupled problem, the implicit relation between pressure and displacement (via static equilibrium) is used to eliminate the displacement field from the pertinent error equations.
Throughout this article, regular font is used to denote scalars (eg, ), bold italic font is used to denote first-and second-order tensors (eg, u, ), and bold font to denote fourth-order tensors (eg, E). The scalar product (single contraction) is denoted by ⋅, double contraction is denoted with : and the outer product denoted by ⊗. For first-order tensors a, b, second-order tensor A and fourth-order tensor B, we thus have for Cartesian components, where repeated indices are summed over (Einstein summation convention). A superposed dot is used for time derivatives (eg,u = du dt ). Volume averaging of a field • is denoted as where Ω □ is the domain occupied by an RVE, and |Ω □ | the corresponding volume. This article is outlined as follows: Section 2 introduces computational homogenization, for the model problem of porous media. Section 3 introduces NMR and describes how it is applied to the microscale (RVE) problem(s). Section 4 discusses how the error in the reduced solution can be estimated in terms of an energy norm, and Section 5 describes the goal oriented error estimation. Section 6 presents numerical results for two example problems, which verify the error estimates. This article is concluded in Section 7, which summarizes the findings.

The model problem-Strong and space-time weak formats
As a model problem, we consider a continuum mechanics description of a linear-elastic porous medium, where the pores are filled with a viscous fluid. We base our model on Jänicke et al 12 which adapts Biot's equations for linear consolidation, 30,31 with displacement u = u(x, t) and pressure p = p(x, t) as the primary fields where is the linear stress tensor, Φ the fluid storage function and w the seepage velocity. The two fields are subjected to standard boundary conditions on the Dirichlet (Γ For simplicity, we will consider linear constitutive relations for the stress, fluid storage, and fluid flux. The stress is given by where E is the constant elastic stiffness tensor, [u] = [u⊗ ] s is the linear strain tensor and is the so-called Biot's coefficient. The storage function and seepage velocity for the liquid phase are given by where is the (initial) porosity, K = kI is the permeability tensor with isotropic permeability k, and is the compressibility parameter of the fluid-filled pore space. and are defined in terms of the bulk moduli of the fluid, K f , and the solid, K s , phase as follows: Finally, we need an initial condition for Φ, viz.
The standard weak format in space-time, corresponding to (3), reads: Find (u, p) ∈  ×  such that The semidiscrete version of (9) is obtained by defining where U h ∈ [H 1 (Ω)] 3 and P h ∈ H 1 (Ω) represents the finite element discretization, and H 1 the space of functions with square integrable derivatives of order 0 and 1. The spaces used in (9) are consequently defined as Here, we introduced Bochner spaces for space-time functions, such that, for example, for a space X of spatial functions with suitable norm || • || X . L 2 (I) and H 1 (I) are the spaces of square integrable functions and functions of square integrable derivatives, respectively, on the time interval I. Remark 1. In order to shorten notation, we henceforth adopt the convention that the trace is included in H 1 (I), meaning that, for example, q ∈ H 1 (I; X) infers that q(0) and q(T) exist and reside in X. □

First-order selective homogenization in the spatial domain
The single-scale problem in Subsection 2.1 is replaced by a two-scale problem. We introduce (i) running averages in the weak form and (ii) scale separation via first-order selective homogenization, following Larsson et al. 32 As a first step, we replace the single-scale problem (9) with the following, two-scale, problem: where the space-time variational forms are defined as In (13), we introduced pertinent space-variational RVE-forms, representing running averages on domains Ω □ located at each macroscale spatial point x: Remark 2. For later use, we also introduce two norms based on □ and (p) □ , respectively: □ In order to define the two-scale trial spaces  FE 2 ,  FE 2 and the corresponding test spaces  FE 2 ,  FE 2 , we introduce first-order homogenization. However, to simplify matters we adopt selective homogenization, in the sense that it is only the displacement field u that is decomposed into macroscale and microscale parts, whereas p is represented only as a fluctuation field. Hence, in each RVE, centered at macroscale coordinate x, we thus decompose u as where we use first-order homogenization for the macroscale part, viz.
with ∶= [u ⊗ ] s . As to the pressure field, we have that is, the pressure field "lives" entirely on the subscale. We are now in the position to define the two-scale ansatz and test spaces cf. Ekre et al. 24 The details of the ansatz and test spaces in (19) is discussed in the following two sections (2.3 and 2.4).

The macroscale (homogenized) problem
The macroscale problem is obtained from Equation (12) by setting v = 0 and q = 0. Find u ∈  such that where  and  are the ansatz and test spaces for the macroscale problem. We omit the exact definitions of these spaces, since we henceforth in this article focus solely on the local microscale RVE-problem. The homogenized stress is defined as where { } is implicit due to the history dependence.

The microscale (RVE) problem
Since we shall (in this article) be concerned only with the solution of the RVE-problem, we consider the situation where u M from (17) is known, i.e. u(t) and (t) are known functions in time (for the given RVE in question). The problem (12) thus reduces to that of finding (u , p) ∈  □ ×  □ that solve where we introduced the RVE space-time variational forms We adopt Dirichlet boundary conditions, for both u and p, and the spaces of spatial functions for the RVE problem are consequently defined as where U □,h and P □,h are the (spatially) FE-discretized function spaces. The trial and test spaces in (22) can be expressed as

Preliminaries
We will now use the fact that Equation (22a) is time-invariant, and does not contain any time derivatives, of u or p, to reduce the problem further. As a preliminary step, we introduce an implicit reduction of the displacement fluctuation u .
For any t ∈ I we define where and the implicit function u p is defined from In summary, for a given RVE, we use the following decomposition of u Remark 3. We note, from (27), that u can be directly assessed in terms of unit strain perturbations, that is, where the "unit fields" □ With the displacement fluctuation implicitly known from p we can formulate a condensed version of the original problem in (22); Find p ∈  □ such that where we defined

NMR-ansatz
We now introduce the concept of NMR with the aim to reduce the computational cost of solving each RVE problem (32).
To this end, we construct a reduced spatial basis for the pressure, and define p R (x,t) as the approximation of p(x,t) with an expansion using N R modes where a=1 is a set of linearly independent basis functions that span the reduced RVE space and where a are mode "activity coefficients." The identification of the spatial modes p a will be discussed further in Subsection 3.4.
Remark 4. We note that in the ansatz (34) the aim is to precompute the spatial modes p a in the "offline" stage, and later solve for the mode activity coefficients in the "online" stage. This is in contrast to other methods, for example, Proper Generalized Decomposition, where both functions are computed in the online stage. □ In order to satisfy (28), and find u p {p(t)}, we make the following ansatz where u a (x) are (spatial) displacement modes and where a (t) are the same mode activity coefficients as those used for the reduced pressure field. The sequence of displacement modes are readily solved from (28). Find u a ∈ U □ such that i.e. one stationary, linear, problem to solve for each mode p a .

Explicit form of the reduced subscale problem
With the approximations from the previous section we can, following the procedure in Jänicke et al, 12 define the reduced equivalent of (22): Find p R ∈  □,R such that where  □,R follows from (25d), with P □ replaced by P □,R . Hence, we can expand the test function q using the spatial pressure modes, i.e. q R = ∑ N R a=1 p a a , and express (38) explicitly as the problem of finding the mode coefficients a (t) ∈ H 1 (I), a = 1, 2, . ., N R such that where Remark 5. Equation (39) can be formulated as a semidiscrete system of size N R as follows where The full two-scale model involves the solution of (reduced) microscale problems (39), typically one for each macroscale quadrature point, in order to find homogenized stress and the sensitivity w.r.t. the macroscale strain . In particular, the homogenized stress can be formulated as where the fourth-order tensor E and the second-order tensors E p,a can be computed in the "offline stage," resulting in efficient evaluation of the homogenized stress, i.e.
where E is the fourth-order stiffness tensor for linear elasticity.

Construction of the reduced basis
In the previous section, we set aside the discussion about the construction of the spatial modes p a for the reduced basis. We merely mentioned that the modes span the reduced pressure space P □,R ⊂ P □ . In this section, we will briefly discuss two different methods to construct the reduced base: POD and SD.

Proper orthogonal decomposition
In order to construct a POD basis for the pressure, we collect N S pressure "snapshots"p k ∶= p(t k ), k = 1, 2, . ., N S , in time, from the solution of the fully resolved problem given in Equation (22). To extract modes from the dataset of snapshots, we first form the correlation matrix C with entries which we use to solve the eigenvalue problem C ⃗ = ⃗, and obtain eigenvalues a and eigenvectors ⃗ a , ordered based on the eigenvalues in decreasing order such that 1 ≥ 2 , … . The final modes, used for the reduced pressure solution, are finally computed as Remark 6. The number of modes N POD R is determined by truncation of the eigenvalues series, in practice when the quotient N ∕ 1 becomes "sufficiently small." □

Spectral decomposition
SD is a well suited method of finding a reduced basis for linear problems. In the context of computational homogenization, a spectral base was used for the transient heat flow problem in Aggestam et al 23 for the RVE problem, and later extended to the full two-scale model in Ekre et al. 24 The problem discussed in this article is not well suited for SD, due to the coupling between the two fields. Nevertheless, it is possible to construct a generalized eigenvalue problem from the uncoupled version by simply ignoring the term where 1 ≤ 2 ≤ … are the eigenvalues, and p 1 ,p 2 , … the corresponding eigenfunctions. From (47), we also obtain the standard orthogonality condition The spectral base is not expected to perform particularly well, since it originates from a modified problem which does not take the coupling into account. However, the spectral base is advantageous for the error estimation, as will be discussed in Sections 4 and 5 and may serve as a "reference solution."

Projection onto the reduced space
For the subsequent error analysis, we define a projection operator Π R ∶ L 2 (Ω □ ) → P □,R such that, for any given q ∈ L 2 (Ω □ ) and for any given reduced set P □,R , the projection Π R q ∈ P □,R is defined by the identity We also define the complementary operator where I is the identity operator. We note that for p, q ∈ L 2 (Ω □ ) the following identity holds

Preliminaries
It is clear that the reduced solution p R is an approximation of p, and that the reduction introduces an error. In this section, we will assess the accuracy of the (reduced) solution by estimating the error. The goal is to obtain an estimate for the error in terms of an energy norm, and in terms of user-defined quantities (goal oriented estimation). We recognize that the total error have multiple sources, in particular (i) NMR and (ii) space and time discretization. In this article, we focus solely on the error from NMR, and thus assume that the resolution in time and space is sufficiently good. The (exact) error in the pressure is defined as where p ∈  □ is the exact solution, and where p R ∈  □,R is the reduced solution. In order to find an estimate for g we use the following "building blocks": • This strategy follows the procedure outlined in Aggestam et al, 23 Ekre et al. 24 This article can be considered a generalization of the previous work, and, in particular, includes the following novelties: the model problem is a coupled two-field problem; the estimator (and the reduced model) are derived for a POD basis (in contrast to a SD basis); the implicit relation between u and p given by static equilibrium is utilized to (i) reduce the error equation to a one-field problem, and (ii) to construct the symmetric norm-inducing form.

The error equation and corresponding residual
From linearity of A □ (•, •) we may establish the error equation: Find g ∈  □ such that where the residual is defined by In the last step, we collected all terms by defining M t and M 0 from the following identities Remark 7. It is important to note that it is not necessary to solve M t explicitly in each time step, since it is enough to compute the sensitivities w.r.t the time-dependent functions. The total computational effort for the error estimate is discussed further in Section 4.5. □ Since we are only concerned with the error stemming from the reduced basis, and hence consider (38) to be solved exactly, we note that and hence obtain the following, Galerkin-like, identity where we in the last step used the projection property from Equation (51).

Auxiliary error equation and associated norm
In order to obtain a norm, and to construct an auxiliary problem where we avoid computation of the implicit relation u p {q} for any q ∈  □ , we introduce the auxiliary bilinear form Â □ (•, •), defined as follows cf. a similar approach in Parés et al. 27 We note that Â □ (•, •) defines a scalar product on the set and hence it can be used to define a norm ||q|| ∶= Furthermore, we note that 3  □ ⊂ □ , and that any functional on  □ can be evaluated for functions in □ . Finally, for any function q ∈  □ , we derive the following inequality i.e. the auxiliary form bounds A □ (•, •) from below. We may now express an auxiliary error equation based on the form in Equation (59). Findĝ ∈ □ such that The norm of the auxiliary errorĝ provides an upper bound of the norm of the true error g which follows from the inequality in (62), the definitions of the error equations in (54) and (63), and from using Cauchy-Schwartz inequality, i.e.

Explicit NMR-error estimates
In what follows, we aim to derive an explicit upper estimate for ||ĝ||, which will thus also serve as an upper bound of ||g||.
From the definition of the norm in (61), we obtain We first note that, from the definition of the auxiliary form Â □ (•, •), the auxiliary error equation (63)

Auxiliary error equation at time t = 0
At time t = 0 the auxiliary error equation (63) reduces to the problem of findingĝ(0) =∶ĝ 0 ∈ P □ such that With Cauchy-Schwarz inequality we obtain the following upper bound from (67): and, thus, obtain an estimate for the contribution to ||ĝ|| at t = 0

Auxiliary error equation at time t ∈ I
In the time interval the auxiliary error equation (63) reduces to the problem of findingĝ(t) =∶ĝ t ∈ P □ such that and we obtain the following upper bound using Cauchy-Schwartz inequality To proceed we need a relation between ||ĝ t || and ||ĝ t || . To this end, we use SD and writeĝ t = ∑ N a=1 a a where a are coefficients to the corresponding eigenmodes a from the solution of the eigenvalue problem with eigenvalues a sorted in increasing order, such that 1 ≤ 2 , … . Withĝ t written as the spectral expansion we can obtain the following inequality Finally, with (71) and (73), we arrive at the following upper bound for the time interval contribution to ||ĝ|| Remark 8. We note that the eigenvalue problem in (72) is equivalent to the eigenvalue problem from Section 3.4.2, where it was used to identify a spectral base for the reduced field. Thus, if such a base is used for the reduction, we may use the same strategy as in Aggestam et al 23 and exchange min to N R , resulting in a sharper bound. □

Auxiliary error equation at time t = T
At time t = T the auxiliary error equation reduces to the problem of findingq(T) =∶q T ∈ P □ such that and, since □ (•, •) is coercive, this leads to the trivial result

Final estimate of the energy norm
Equation (66) The estimate is similar to the estimate given for the transient heat flow problem in Equation (112) from Aggestam et al, 23 with some important differences: • in Aggestam et al 23 the norm inducing functional is equal to the original functional for the same argument 4 , whereas here, in Equation (65), we need an extra inequality, ||g|| 2 = Â □ (g, g) ≤ A □ (g, g), due to the handling of the term coupling the fields; • in Aggestam et al 23 the eigenvalue in the estimate is chosen to the largest resolved eigenvalue, which is possible since a spectral based is used for the reduction, whereas here we need to use the lowest eigenvalue in the general case. 5

Quantification of computational effort for the estimate
The final estimate in Equation (77)  In particular it is important to note that the number of auxiliary solves is fixed and does not scale with neither mesh size nor the number of time steps. These 6 + N R functions can be precomputed, and, when they are known, the estimate in (77) is fully explicit and efficient to evaluate.

Preliminaries
In this section, we will discuss quantification of the error in terms of user-defined quantities of interest (QoI) using goal oriented error estimation. The procedure to obtain the explicit estimate follows the steps in Section 4 closely, and in particular, we will again only focus on the error introduced by the reduced base, and ignore space-time discretization errors by considering the fully resolved finite element problem to be the exact solution.
For the goal oriented error estimation, we use the following "building blocks": • definition of linear output functional for the chosen quantity of interest (Section 5.2); • definition of the dual problem and corresponding reduced dual problem (Section 5.3); • formulation of the dual residual and the dual error with upper and lower bounds (Sections 5.4 and 5.5); • construction of explicit error bounds for the error in the quantity of interest (Section 5.6).

Linear output functional
We introduce linear goal functionals Q (u) □ , corresponding to the original (uncondensed) RVE problem in (22), that represents the Quantity of Interest (QoI) □,T (p(T)).
(78b) 4 In particular, from Equation (74) in Aggestam et al 23 there is ||e|| 2 = A s □ (e, e) = A □ (e, e) for the error e. 5 As mentioned in Section 3.4.2 it is possible to use a spectral base for problem discussed in this article, in which case the largest resolved eigenvalue can be used in the estimate, as noted in Section 4.4.2. 6 The auxiliary problems are of size N, for example, the size of the underlying finite element problem that the reduced basis is based on, and have a "mass matrix" as the left-hand side.
For the condensed problem (32), we can clearly reformulate the QoI as However, for the subsequent error analysis we wish to establish u p {•} for the entire P □ , and not just for P □,R as was done in Subsection 3.2. To this end, we establish the following problems to solve for auxiliary "influence functions" u ⋆ t (t) ∈ U 0 □ and u ⋆ T ∈ U 0 □ : We may now rewrite (79) as an explicit expression in p, by using the auxiliary dual problems in (80), and the definition of the implicit function u p in (28) Equation (81) now represents the output and is fully explicit in p. Finally, we may express the functional for the error which is an explicit functional on the space

Example: Time-averaged homogenized stress
Consider the ij-component of the homogenized stress, ij , as the Quantity of Interest. The output functional can thus be defined by We note, in particular, that this functional does not depend on time, and therefore it is enough to solve the auxiliary problems (80a) once, instead of once per time step, which is needed in the general case.

Dual problem
We define the dual problem as that of finding p ⋆ ∈  □ such that where we recall the test space  ⋆ □ from (83), and where the dual form A ⋆ □ is defined as We note that, for q, r ∈  □ which follows from integration by parts, namely, and The reduced dual problem is defined as that of finding p ⋆ R ∈  □,R such that where  ⋆ □,R follows from (83), with P □ replaced with P □,R .

Dual error equation and corresponding residual
We define the dual error equation as the problem of finding g ⋆ ∈  □ such that where g ⋆ ∶= p ⋆ − p ⋆ R is the error in the dual solution. Consequently, due to (87) and the Galerkin orthogonality in (54) we obtain since p ⋆ R ∈  □,R ⊂  □,R . The dual residual is expanded similar to Subsection 4.2 where we introduced Consider now the following equalities 7 : where we in the last step used the following identities, derived from the error equations and the Galerkin orthogonality In a procedure resembling the Parallelogram law, cf. Oden and Prudhomme, 29 we subtract the "−"-equation from the "+"-equation in (104) and obtain where the terms in the right-hand side of (106) can be bounded from below, and from above using Cauchy-Schwartz inequality, We can now obtain upper and lower bounds on Q □ (g) by combining (106) with (107) and (108) In what follows we will use the trivial lower bound, ||g ± || ≥ 0. Utilizing a sharper lower bound for ||g ± || would naturally result in a sharper estimate. However, it requires implicit procedures, which would make the final estimate more computational demanding. As a last step, we also utilize the bounds given by (65) and (98) and arrive at upper and lower bounds for the error in the quantity of interest

Explicit bounds for the error in quantity of interest
In order to obtain explicit upper and lower bounds for E − Q and E + Q from Equation (110) we need an estimate for the norm of the composite error, ||ĝ ± ||. To this end, we use the strategy described in Subsection 4.4 and consider the composite error equation (101) localized in time for times t = 0, t ∈ I and t = T separately.

5.6.1
Composite error equation at t = 0 At time t = 0 the composite error equation reduces to findingĝ ± (0) =∶ĝ ± 0 ∈ P □ such that In terms of the norm, we obtain the following upper bound, using Cauchy-Schwartz inequality and obtain the following contribution to ||g ± || at t = 0

Composite error equation for t ∈ I
In the time interval, t ∈ I, the composite error equation reduces to findingĝ ± (t) =∶ĝ ± t ∈ P □ such that and, by using the same arguments as in Section 4.4.2, we obtain the following upper bound for the contribution to ||ĝ ± ||

Composite error equation for t = T
At time t = T the composite error equation reduces to findingĝ ± (T) =∶ĝ ± T ∈ P □ such that and with the same arguments as in Section 5.6.1 we obtain the following upper bound of the contribution to ||ĝ ± ||

Final upper and lower estimate for the error in quantity of interest
With the norm from Equation (66), the bounds from Equations (113), (115), and (117), and the upper and lower bounds of the Quantity of Interest from Equation (110) we obtain the following explicit estimates Finally, the estimates can be optimized by finding an optimal value of . We note that Equation (118) can be written as Solving the trivial minimization problem, we obtain 2 opt = √ C A , which results in the following explicit upper and lower bounds Once again we emphasize that these are all computations that can be performed in the "offline" stage, stored, and then used for the error estimation throughout the "online" stage.

NUMERICAL EXAMPLES
For the numerical examples we consider two different RVEs in three spatial dimensions. The RVEs consist of gas saturated matrix material (phase 1) with spherical, water saturated, inclusions of two different types (phases 2 and 3, respectively). The material parameters for the three phases are presented in Table 1.
In order to quantify the sharpness of the error estimate (as compared with the exact error) we define the effectivity index as where E est is the estimate (Equation (77)), and E the exact error. Similarly, for the goal oriented estimate we define where E Q,est is a "worst case" combination of the upper and lower bounds, and E Q the exact error in the quantity of interest.

Example 1: Isotropic RVE with a single centered inclusion
As a first example, we consider an isotropic RVE of size 1 m × 1 m × 1 m, with a single centered inclusion with radius 0.4 m, see Figure 1. For this example, we use the same load case as in the training computations (see Section 6.1.1) for the prediction; the 11 component is rapidly increased and then kept constant, while keeping the other strain components at 0.

Example 1: Identification of reduced basis
The snapshots for the POD are extracted from training computations. For the first example, we perform a stress-relaxation test as training by ramping up 11 (t) = (t) according to Figure 2, while keeping the other components at 0. The resulting homogenized stress is plotted in Figure 3. The snapshots for the POD are collected from the "relaxation part", i.e. from the time steps where t > 1 × 10 −5 s. In Figure 4 the first two POD modes are plotted, together with the first SD modes for the RVE.

Example 1: Energy norm of the error
The estimated error in energy norm, E est , given in Equation (77), is plotted in Figure 5 for the prediction load case (which, for this example is the same as the training load case). The exact error, E = ||g||, is also plotted in the same figure. We note that, as expected, the residual, and therefore also the error estimate, decreases much faster for the POD base compared with the SD base. The effectivity index, , for the estimate in energy norm is plotted in Figure 6. We note that the effectivity index increases quickly for the POD base, but stays rather constant for the SD base. One explanation for this is the fact that, with the SD base, we can use a higher eigenvalue for the final estimate, as discussed in Section 4.4.2.

Example 1: Error in time-averaged homogenized stress
In Section 5, we discussed how to estimate the error in a user-defined quantity of interest, and in particular the homogenized stress. We now consider the time-averaged homogenized stress, 11 , as the quantity of interest and compute the upper and lower estimates for the error, E + Q,est and E − Q,est , given in Equation (121). The results are plotted in Figure 7, and the combined estimate, from Equation (123), is plotted in Figure 8. In both figures, we also include the exact error, E Q , as compared with the full finite element solution. We note that the behavior of the estimates, and the corresponding effectivity indexes in Figure 9, is very similar to the estimate of the energy norm from the previous section. However, in Figure 9, we see that the effectivity index deteriorates even more for the POD estimate.

Example 2: Heterogeneous RVE with multiple inclusions
For the second example, we consider an RVE of size 1 m × 1 m × 1 m, with multiple, nonoverlapping, spherical inclusions with radii between 0.1 and 0.4 m, see Figure 1. For the prediction, we combine prescribed macroscale strain in 11, 22, and 33 directions according to Figure 10, resulting in macroscale stresses presented in Figure 11.

Example 2: Identification of reduced basis
For the second example, we also use the stress relaxation test from Section 6.1.1, but now in three directions by setting ii (t) = (t), for i = 1,2,3, resulting in three training computations from which we collect snapshots. The resulting first two POD modes are plotted in Figure 12, together with the first two SD modes.

Example 2: Energy norm of the error
The estimated error in energy norm, E est , given in Equation (77), is plotted in Figure 13, together with the exact error, E = ||g||. We note that the result is similar to the results for the first example in the previous section. In particular, we note that, as expected, the residual, and therefore also the error estimate, decreases much faster for the POD base compared with the SD base. The effectivity index, , for the estimate in energy norm is plotted in Figure 14. We note that the effectivity index is higher for the POD base than for the SD base, just as for the previous example.

Example 2: Error in time-averaged homogenized stress
Once again we use the time average of 11 as the quantity of interest for the goal oriented error estimation. The resulting upper and lower bounds, E + Q,est and E − Q,est from Equation (121), are plotted in Figure 15, together with the combined estimate and effectivity index (123) plotted in Figures 16 and 17, respectively. The results are similar to the first example, where the POD base results in much smaller residual and hence smaller estimate, but the spectral base results in a better effectivity index.

Example 3: Influence of training data
In the third numerical example, we investigate the effect of the training data on the error estimate. The training data are collected using the same method as for the previous example (see Section 6.2.1), however, we vary the activation function F I G U R E 18 Example 3: Threshold number of modes N R required for reaching 10%, 1%, 0.1%, and 0.01% estimated relative error in terms of the energy norm, respectively, for the three different POD basis constructed from the datasets DS1, DS2, and DS3 in order to obtain three different datasets (DS1, DS2, DS3) to build the POD basis from. The microstructure of the RVE, and the load case used for the prediction, are the same as in Example 2.
In Figure 18, we illustrate the threshold number of modes needed in order to obtain a certain accuracy, in terms of the energy norm, indicated by the estimate. We note that the threshold differs between the different datasets quite a lot even for this simple load case.
This example also illustrates the need for different number of modes, or even different basis, for different macroscopic location in an FE 2 setting where the loading are different for each RVE computation.

CONCLUSIONS AND OUTLOOK
In this article, we have presented an a posteriori error estimate for the error introduced by NMR in the context of homogenization of porous media. More specifically, we extended the reduced order model presented in Jänicke et al 12 with an explicit, fully computable, and computationally efficient error estimator for guaranteed bounds of the error in terms of (i) energy norm and (ii) user-defined QoI. The estimator is a generalization of similar estimators presented by the authors 23,24 with the major important differences that it handles the coupled problem, and that the estimator is derived for a POD basis rather than an SD basis. The performance of the estimator was demonstrated by two numerical examples, and for reduced basis obtained by (i) SD and (ii) POD. As expected, the POD base results in much smaller residuals compared with the SD base, but in general the effectivity index of the estimator when using an SD base is better. It is noted that the sharpness of the estimator quickly deteriorates with increasing number of modes. However, considering the fact that the estimates are fully computable with negligible overhead, a sharper estimate would probably require significantly higher computational cost.
An interesting extension, left for future work, would be to develop a method to create a mixed base by combining spectral modes and POD modes in a clever way, in order to reduce the high effectivity index obtained from the estimator for the pure POD case.
As to future work, the obvious next step is to embed the error estimator for the subscale problem, presented in this article, into a full-fledged nested FE 2 procedure. This would require taking care of error transport between the two scales, similar to what was presented in Ekre et al 24 for the linear heat flow problem. We also, once again, note that in this work we completely ignored errors from sources such as time and space discretization, and it would be of interest to extend the estimator with the ability to estimate those errors as well. Finally, we note that we have not considered adaptive techniques in this article. It would be interesting to investigate this, and use the error estimator as a measure for adaptively adding new modes to the subscale problem on the fly.