Combining spectral and POD modes to improve error estimation of numerical model reduction for porous media

Numerical model reduction (NMR) is used to solve the microscale problem that arises from computational homogenization of a model problem of porous media with displacement and pressure as unknown fields. The reduction technique and an associated error estimator for the NMR error have been presented in prior work, where both spectral decomposition (SD) and proper orthogonal decomposition (POD) were used to construct the reduced basis. It was shown that the POD basis performs better w.r.t. minimizing the residual, but the SD basis has some advantageous properties for the estimator. Since it is the estimated error that will govern the error control, the most efficient procedure is the one that results in the lowest error bound. The main contribution of this paper is further development of the previous work with a proposed combined basis constructed using both SD and POD modes together with an adaptive mode selection strategy. The performance of the combined basis is compared to (i) the pure SD basis and (ii) the pure POD basis via numerical examples. The examples show that it is possible to find a combination of SD/POD modes which is improved, i.e. it yields a smaller estimate, compared to the cases of pure SD or pure POD.


Introduction
Multiscale modeling with computational homogenization is a well-known approach for material modeling. The main advantage over direct numerical simulation (DNS) is reduced computational cost. One standard method for multiscale modeling is the finite element squared (FE 2 ) procedure, where subscale computations are carried out on representative volume elements (RVE) in each point of the macroscale domain 1 in a nested iteration scheme. Even though there are benefits with FE 2 over DNS, the FE 2 scheme is still very computationally demanding for practical problems, especially for fine macroscale meshes in three dimensions, where the number of RVE problems rapidly increases with mesh density. It is therefore of interest to reduce the computational cost of solving the individual RVE problems. 1 In practice one RVE solution per macroscale quadrature point in a finite element setting.
B Fredrik Ekre fredrik.ekre@chalmers.se 1 Department of Industrial and Materials Science, Chalmers University of Technology, 41296 Gothenburg, Sweden A number of numerical model reduction 2 (NMR) methods have been proposed for reducing the solution space of a discrete RVE problem. We highlight, in particular, methods based on superposition of "modes" that are characteristic to the solution field. Waseem et al. [1] and Aggestam et al. [2] presented reduced models for computational homogenization of linear transient heat flow based on spectral decomposition (SD). In the context of small strain (visco)plasticity various attempts have been made to approximate the inelastic strains with "inelastic modes". One example is the "eigendeformation-based-reduced-order homogenization" technique introduced by Fish and coworkers [3,4], which relies on transformation field analysis (TFA), 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 visco-elasticity and a class of stan-dard dissipative materials. Jänicke et al. [12,13] applied this approach to poroelasticity, where 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, and the "mode coefficients" can be interpreted as internal variables.
Naturally the use of a reduced basis results in an additional source of error which is of interest to control and quantify. 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. [14,15] for heterogeneous multiscale methods and by Boyaval [16] 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 have been proposed by Kerfriden et al. [17] and Chamoin and Legoll [18]. Aggestam et al. [2] presented an estimator, with guaranteed bounds for the NMR error, for the SD-reduced RVE problem pertaining to transient linear heat flow. The estimator is based on an auxiliary symmetric form of the original problem, cf. Pares et al. [19][20][21] who presented this strategy for estimation of discretization errors, and fully computable bounds on the NMR error were derived based on the discrete residual, cf. e.g. Jakobsson et al. [22]. The estimator was derived for (i) an energy norm and (ii) for arbitrary user-defined quantities of interest in a procedure inspired by Oden and Prudhomme [23]. The estimator presented in [2] was generalized for a POD basis in the context of poroelasticity by Ekre et al. [24] where it was shown that the POD basis outperforms the SD basis when it comes to minimizing the residual (and hence the residual-based estimator), but the SD basis behaves better in that it does not overestimate the error as much.
For the NMR estimator in Aggestam et al. [2] and Ekre et al. [24] the fully resolved finite element solution is considered exact, i.e., the estimator only quantifies the error stemming from the use of a reduced base and "ignores" discretization errors. Naturally an estimator which quantifies these errors as well would be desirable, 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, many of the necessary quantities can be precomputed in an "offline" stage where one can afford a fine mesh without having to tackle the increased cost of integration, as one would have to do for e.g. a non-linear problem.
This paper is a continuation of the work presented in [12,24]. In particular, we consider NMR for a continuum mechanics model of porous media, cf. Jänicke et al. [12], and an associated error estimator, cf. Ekre et al. [24]. The main contribution is a proposed combined basis with the inten-tion of combining the "residual minimizing" property of the POD basis with the "estimator improving" property of the SD basis, and in particular investigate whether it is possible to find an improved composition. In terms of error control, the solution will be accepted as sufficiently accurate when the estimated error is below the prescribed tolerance. Hence, at least when considering guaranteed bounds of the error, the most efficient sequence of approximations is the one that has the best convergence in terms of the estimated bound on the error. In other words, the procedure that generates the lowest (exact) error is not necessarily the most practical one if the estimated error (that has to be used for the pertinent error control) is large.
Throughout this paper, regular font is used to denote scalars (e.g. α), bold italic font is used to denote first and second order tensors (e.g. u, ε), and bold font to denote fourth order tensors (e.g. E). The scalar product (single contraction) is denoted with '·', double contraction is denoted with ':' and the outer product is denoted with '⊗'. 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 (e.g.u = du dt ). Volume averaging of a field • is denoted as where is the domain occupied by an RVE, and | | the corresponding volume.
The remainder of this paper is outlined as follows: Sect. 2 introduces computational homogenization for the model problem of porous media. Section 3 introduces numerical model reduction (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 (i) an energy norm and (ii) user-defined quantities of interest. Section 5 presents numerical results for two example problems, which verify the error estimates, and Sect. 6 summarizes and concludes the paper.
2 Two-scale analysis based on computational homogenization

The model problem: strong format of linear-elastic porous media
As a model problem we consider a continuum mechanics description of a linear poroelastic medium, where the pores are filled with a viscous fluid. We base our model on Jänicke et al. [12] which adapts Biot-Willis equations for linear consolidation [25,26], with displacement u = u(x, t) and pressure p = p(x, t) as the primary fields satisfying where σ is the Cauchy 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 kinematics symmetric 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 = k I is the permeability tensor with isotropic permeability k, and β is the effective 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.

First order selective homogenization in the spatial domain
In order to derive the pertinent two-scale formulation we follow the standard procedure of computational homogenization, see e.g. Larsson et al. [27]. For brevity we summarize the steps here and refer to Ekre et al. [24] for details. To obtain the macro/microscale equations from the strong form we: • define the weak space-time format of Eq. (3); • homogenize the spatial integrals by introducing running averages over RVEs, located at macroscale pointsx; • introduce scale separation and first order selective homogenization 3 , e.g., whereū is the macroscopic displacement field, and where u μ and p μ are the microscopic (fluctuation) displacement and pressure fields, respectively; • derive the macroscale problem by adopting the pertinent macroscale testfunction v M (v μ = 0); • derive the microscale problem by adopting the pertinent microscale testfunction v μ (v M = 0), one for each RVE.

The macroscale (homogenized) problem
The resulting macroscale problem reads; Findū ∈Ū such that whereŪ andV are the ansatz and test spaces for the macroscale problem. We omit the exact definitions of these spaces, since we henceforth in this paper 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
In this paper we are concerned only with the solution to the RVE-problem, and thus consider the situation where u M from (9a) is known, i.e.ū(t) andε(t) are known functions in time (for the given RVE in question). The resulting microscale problem then reads; Find (u μ , p) ∈ U μ × P that solve where we introduced the RVE space-time variational forms and the space variational forms representing the running averages over each RVE Remark For later use we also introduce two norms based on m and a ( p) , respectively: 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. We thus consider the fluctuation of u and the pore pressure itself to vanish on the boundary. The Bochner trial and test spaces 4 in (12) can be expressed as P :=H 1 I ; P , (17c) Q :={q(x, t) : q| t=0 ∈ P , q| I ∈ L 2 (I ; P )}. (17d)

Preliminaries
As a preliminary step we follow [24] and utilize the timeinvariance of Eq. (12a) to introduce an implicit reduction of the displacement fluctuation u μ . For any t ∈ I we define where u μ ε (t) ∈ U 0 , and the implicit function u μ p , are defined such that they fulfill (12a). With the decomposition in (18) we thus consider the following decomposition of u within each RVE: Finally we formulate a condensed version of the original problem in (12); Find p ∈ P such that where we defined Equation (20) is the starting point for the numerical model reduction and the error estimation.

NMR-ansatz
We now introduce numerical model reduction (NMR) with the goal of reducing the computational effort of solving each RVE problem (20). 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). For the approximation we use N R modes as follows where { p a } N R a=1 is a set of linearly independent basis functions that span the reduced RVE space and where ξ a are mode "activity coefficients". In previous work by the authors [24] the spatial modes have been identified using (i) spectral decomposition (SD) and (ii) proper orthogonal decomposition (POD). In this paper we propose to combine SD modes and POD modes in an attempt to minimize the error bound. We therefore define N SD R and N POD R as the number of SD modes and POD modes, respectively, and expand the approximation (22) in terms of these different modes is the set of POD modes, and where ξ SD a and ξ POD a are the corresponding mode activity coefficients. We adopt the following enumeration of the modes and similarly for the mode activity coefficients and henceforth use the compact notation from Eq. (22) when differentiation between SD modes and POD modes is not necessary. SD modes are obtained from Eq. 12b (by ignoring the coupling). POD modes are obtained by performing representative, fully resolved, training simulations to collect pressure snapshots which then are used as the basis for the POD procedure. For the details of the mode identification we refer to Ekre et al. [24].

Remark
We note that all p S D a are linearly independent, and that all p P O D a are linearly independent, by construction. However, this does not guarantee that a combination of basis is linearly independent and, thus, it is important to ensure this property in the combining procedure. Alternatively it would be possible make orthogonalization part of the basis construction procedure using, e.g., Gram-Schmidt orthogonalization.
We utilize linearity of the sensitivity u where each displacement mode is solved from each corresponding pressure mode under the constraints of (12a): 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 (12): Find p R ∈ P ,R such that where Q ,R follows from (17d), 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 (29) explicitly as the problem of finding the mode coefficients ξ a (t) ∈ H 1 (I ), a = 1, 2, . . . , N R . The RVE problem is now reduced to a semi-discrete system of size N R where (31d)

Preliminaries
Using a reduced basis for the solution obviously results in an approximate result, since p R = p. In this section, we will present a method for assessing the accuracy of the reduced solution and approximate the error. We note that the error has other sources besides the NMR. In particular the error is also a result of the time-and space-discretization used for the solution. In this paper we focus entirely on the NMR error and assume that the discretization errors are negligible compared to the NMR error. Effectively this means that we consider the fully resolved finite element solution to be exact, i.e. that a sufficiently fine discretization is used. We thus consider p = p h ≈ p R , where p h is the discrete solution and define the error as The following "building blocks" will be used in order to derive fully explicit and computable error estimates for the reduced basis:

The error equation and corresponding residual
From linearity of A (•, •) we may establish the error equation: Find g ∈ P such that where is the exact error and where the residual is defined by where the Galerkin orthogonality-type identity R (q) = R ( C q) follows from (i) ignoring space-time errors, 5 and (ii) the definition of the projection operator C = I − R .
Here the projection into the reduced set of spatial functions is defined as follows: M t and M 0 are defined by collecting all the right hand side terms in m (•, •) using the following identities Remark 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.

Linear output functional, dual problem and dual error equation
In order to measure the error in arbitrary, user-defined, quantities of interest (QoI) we introduce linear goal functionals Q (u) , Q ( p) , corresponding to the original (uncondensed) RVE problem in (12) For the condensed problem (20), we can clearly re-formulate 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 Sect. 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 (38) as an explicit expression in p, by using the auxiliary dual problems in (39), and the definition of the implicit function u p (see Sect. 3.1) ,T ( p(T )) .
(40) Equation (40) now represents the output and is fully explicit in p. Finally, noting that Q ( p) is affine, we may express the linear functional ,T (q(T )), , Example 2 Homogenized stress at t = T Similarly to the example above, the homogenized stress at t = T as a quantities of interest can be obtained using the following definitions We note, in particular, that these example functionals does not depend on time, and therefore it is enough to solve the auxiliary problems (39a) once, instead of once per time step, which is needed in the general case.
We now define the dual problem and the corresponding reduced dual problem as that of finding p ∈ P such that where we recall the test space Q from (42), from which we also obtain Q ,R by replacing P with P ,R . The dual form A is defined as We can now define the dual error equation: Find g ∈ P such that where g := p − p R is the error in the dual solution and where the dual residual is given by where the identity R (q) = R ( C q) follows from the same argumentation as in Sect. 4.2 i.e. the solution to Eq. (45b) is considered to be exact. M t and M T are defined by collecting all terms in m (•, •), i.e. (49b) Remark We once again note that it is not necessary to solve for M t explicitly in each time step, it is enough to solve for the sensitivities w.r.t. the time dependent functions.

Auxiliary bilinear form, error equations, and associated norm
Following the strategy in Ekre et al. [24] we introduce an auxiliary bilinear form cf. the work by Parés et al. [21]. SinceÂ (•, •) defines a scalar product on the spacê P :={q(x, t) : q| t=0 ∈ P , q| I ∈ L 2 (I ; P ), it induces a norm that will be used to estimate the error where we also note the bound w.r.t. the symmetric part of the original bilinear form.
With the definition of the auxiliary form we also define the following two error equations for auxiliary error representationsĝ andĝ Following from the primary and dual error equations, in (33) and (47), the corresponding auxiliary error equations in (53), the relation between A andÂ in (52), and the Cauchy-Schwartz inequality, we find that the norm of the true errors g and g are upper bounded by the norm of the corresponding auxiliary error representationsĝ andĝ , viz. g ≤ ĝ , and g ≤ ĝ . (54)

Preliminaries
In order to obtain explicit and fully computable estimates in terms of (i) the energy norm and (ii) arbitrary quantites of interest expressed in the linear goal functional we will first consider the following generic error equation; Find χ ∈P s.t.
where the residual is defined as for generic functions M t (x, t) ∈ P and M 0 (x), M T (x) ∈ P . Given this configuration we derive an explicit and fully computable upper estimate of χ , see "Appendix A", given by where λ N SD R is the largest spectral eigenvalue used for the reduced basis.

NMR error in energy norm
In order to obtain an estimate for the NMR error in terms of an energy norm we apply the generic error estimate from (57), and choose χ =ĝ, resulting in in the residual. The estimate for the energy norm of the error thus becomes

NMR error in quantity of interest
In order to find explicit NMR-error estimates for the quantity of interest we first define the "composite errors", g ± andĝ ± , as linear combinations of the corresponding primal and dual parts, viz.
g ± := κ g ± κ −1 g ∈ P ,ĝ ± := κĝ ± κ −1ĝ ∈P , for any κ = 0. From linearity, it follows that 6 A (g ± , δg) = R ± (δg) ∀δg ∈ P , and where Consider now the following equalities, which follow from the error equations and the Galerkin orthogonality, cf. [24], Similarly to the procedure in Oden and Prudhomme [23] we substract the "−"-equation from the "+"-equation and obtain To obtain upper and lower estimates we note that the terms on the right hand side can be bounded from below, using Eqs. (33), (47), (52), and (61), and from above using the Cauchy-Schwartz inequality, By combining Eq. (63) with Eqs. (64), (65), the trivial bound g ± ≥ 0, and the relation between • and • from Eq. (54), we obtain the following upper and lower bound on the quantity of interest: The last step is to apply the generic estimate from Eq. (57) to obtain estimates for ĝ − 2 and ĝ + 2 , respectively. We note that the auxiliary composite error can be solved from the generic error equation in (55), with χ =ĝ ± , and with The final upper and lower bounds on the error in quantity of interest becomē 6 Note that we here use the subscript ± as an alternative to presenting two equations; one "+" and one "−".Q (68b)

Adaptive mode selection strategy
For the purpose of obtaining the lowest estimate it is of interest to find a strategy for adaptive mode selection, i.e. a strategy for determining if the next mode should come from the pool of SD modes or the pool of POD modes. The modes in each pool are ordered based on relevant eigenvalues. Hence, we only have to determine which out of the two candidates to include in the approximation space. For a basis with N R modes with estimate E N R est we want to predict E N R +1 est for the two alternatives. We note that, ignoring the boundary term, the estimate in (59) can be written as The results from Ekre et al. [24] suggest that the strength of the SD modes is mainly that a higher eigenvalue can be used in the estimate. In order the estimate E N R +1 est , for the case of one extra spectral mode, we assume that its contribution is only the increased eigenvalue, and ignore any effect this mode have on the residual, and obtain i.e. with one extra SD mode the estimate is approximated to reduce with the factor φ SD .
To estimate E N R +1 est , for the case of one extra POD mode, we assume that the estimate is reduced with the same factor as the previous POD mode addition, i.e.
i.e. with one extra POD mode the estimate is approximated to reduce with the factor φ POD . We can now compare φ SD and φ POD and choose the next mode based on the estimated error estimate reduction. Finally, in order to not get stuck on local "plateaus" (e.g. consecutive eigenvalues that are close in the spectrum, or "bad" POD modes), if both φ SD and φ POD are large, e.g. φ SD > ρ and φ POD > ρ for some value ρ, we instead alternate between SD and POD modes. Remark Clearly, at the extreme case of using all available POD modes, enrichment with SD modes remains as the only option.

Numerical examples
In this Section we present numerical results that showcase the estimator. As a performance measure of the estimator, as compared to the exact error, we define the effectivity indices η and η Q for the estimator of the energy norm and quantity of interest, respectively where E est and E are the estimate and exact error in energy norm, and where E Q,est and E Q are the ("worst case") estimate and exact error in the quantity of interest. The exact errors E and E Q are obtained by computing reference finite element solutions using the original discretization that was used for obtaining the reduced basis.

RVE configurations and eigenvalue spectra
We consider three different RVE configurations in three spatial dimensions, see Fig. 1. The RVEs consist of gas saturated matrix material (phase 1) with spherical, water saturated, inclusions with various permeability (phases 2 and 3, respectively). The material parameters are presented in Table 1.
The experimental setup is the same as in Ekre et al. [24] and training computations for the POD modes are performed in the same manner, e.g. by simulating stress relaxation after a rapid increase of strain in the 11, 22, and 33 directions, respectively, cf. Sections 6.1.1 and 6.2.1 in [24]. We note from the previous Section that, in order for the estimator to be improved by spectral modes, the eigenvalue sequence is important, and in particular for the first eigenvalues since the region of interest is where N R N . The first 25 eigenvalues for the three RVEs are plotted in Fig. 2. We note that the eigenvalues for RVE-1 and RVE-3 increase much faster than for RVE-2; for RVE-1 a factor of 10, compared to the first eigenvalue, is reached already at the second eigenvalue, whereas for RVE-2 a factor 10 is reached only after 20+ modes. Based on this we can predict that the addition of spectral modes to the POD basis will be more beneficial for RVE-1 than for RVE-2. The spectrum for RVE-3 is similar to that of RVE-1 and both of them reach a factor 100 around 25 modes, however, RVE-3 has multiple consecutive identical values meaning that several spectral modes need to be added before the estimator benefits from an increased eigenvalue.

Example 1: Heterogeneous RVE with multiple inclusions
For the first example we consider an RVE of size 1m × 1m × 1m, with multiple, non-overlapping, spherical inclusions with radii between 0.1m and 0.4m, see RVE-1 in Fig. 1.
The RVE is prescribed with macroscale strain history in 11, 22, and 33 directions according to Fig. 3, which results in macroscale stresses presented in Fig. 4.

Energy norm of the error
The estimated error in energy norm, E est (Eq. (59)), and the corresponding exact error, E = g , are plotted in Fig. 5 for multiple solutions with varying number of spectral modes included in the reduced basis. The corresponding effectivity index, η, is plotted in Fig. 6. We note that, as the combined basis includes more spectral modes, the estimate performs better 7 . In Ekre et al. [24] it was shown that POD modes are superior, compared to SD modes, in terms of minimizing the residual, and, thus, the improvement for the combined basis is attributed (almost entirely) to the larger eigenvalue used in the estimate. At a certain threshold, however, it will be more benefitial to enhance the basis with only POD modes since the eigenvalue spectrum typically flattens out as seen in Sect. 5.1.
In order to understand the advantage of the combined estimator it is helpful to study the behavior for a fixed total number of modes, with varying composition of SD and POD modes. In Figs. 7 and 8 the estimate, exact error, and corresponding effectivity index is plotted for N R = N SD R + N SD R = 5 and N R = N SD R + N SD R = 10, respectively. From the figures it is clear that, for a given number of total modes, there is an optimal number of spectral modes. Note also that, although the effectivity index descreases monotonically with number of spectral modes, the raw value of the estimate starts to increase after a certain threshold number of spectral modes.
In Fig. 9 the result of the adaptive procedure (cf. Sect. 4.6 with ρ = 0.9) is shown. For this example the adaptive procedure results in a basis that give a smaller estimated error compared to the pure SD or pure POD basis.

Error in time-averaged homogenized stress
In Sect. 4 we discussed how to estimate the error in a user-defined quantity of interest, and in particular the homogenized stress. For this example we consider time-averaged homogenized stress,σ 11 , as the quantity of interest and compute the resulting combined estimate, E Q,est from Eq. (72b). The result is plotted in Fig. 10 together with the corresponding exact error, E Q . The behavior of the estimates, and the corresponding effectivity indices in Fig. 11, are similar to the estimate of the energy norm presented in the the previous section-including spectral modes is advantageous, up to a certain threshold, for the estimator, and result in a better effectivity index.
In Figs. 12 and 13 the estimated and exact error is plotted for a fixed total number of modes; N R = N SD R + N SD R = 5 and N R = N SD R + N SD R = 10, respectively. Similar to the energy norm estimate we note that for a given total number of modes there is an optimal composition of SD and POD modes that minimizes the estimate. 7 Note that for N SD R = 0, i.e. a pure POD basis, we still need to compute the first eigenvalue for the estimator. Thus, N SD R = 0 and N SD R = 1 use the same eigenvalue for the estimate.

Example 2: RVE with a single centered inclusion
As a second example we consider an RVE of size 1m × 1m × 1m, with a single centered inclusion with radius 0.4m, see Fig. 1. From Fig. 2 it is evident that the eigenvalue spectrum for this RVE is much less advantageous for the estimator since the initial eigenvalues are very close. In addition, for this example we use the same load case as in the training computations, cf. Ekre et al. [24], and simulate stress relaxation after a rapid increase of normal strain.

Error in time-averaged homogenized stress
Once again we consider time-averaged homogenized stress, σ 11 , as the quantity of interest and compute estimate, E Q,est (Eq. (72b)), and exact error E Q . The results are plotted in Fig. 15. For this example, the advantage with the higher eigenvalue in the combined basis is almost entirely diminished and using additional SD modes is comparable to using a pure POD basis. Note, however, that the effectivity index, plotted in Fig. 16 is still noticeably improved. Figure 14 shows the result of the adaptive procedure. Similar to the energy norm estimate from the previous section, the combined basis give a smaller estimated error compared to the two cases of pure SD or pure POD basis.
Once again we plot the estimate and the exact error for a fixed total number of modes, in Fig. 17 for N R = N SD R + N SD R = 5 and in Fig. 18 for N R = N SD R + N SD R = 10. Compared to the previous example, where there were a distinct minimum for the estimator for N SD R > 1, the benefit of enhancing the basis with SD modes, rather than POD modes, is almost entirely diminished since the POD modes decrease the residual with a higher rate than the eigenvalues increase. Figure 19 shows the result of the adaptive procedure. For this example the first SD mode is not added until mode number eight when the effectiveness of additional POD modes decrease.

Example 3: RVE with an array of spherical inclusions
In the third example we consider an RVE with multiple spherical inclusions, with varying permeability k, arranged in a 3 × 3 × 3 grid, see Fig. 1. The eigenvalue spectrum for this configuration is plotted in Fig. 2. It is in particular of interest to note that the first four eigenvalues are very close to each other, and it is thus necessary to add five spectral modes to the basis before the estimator benefits from an increased eigenvalue. The training and loading is the same as for Example 1.

Error in homogenized stress at time t = T
We now consider homogenized stress at time t = T ,σ 11 (T ), as the quantity of interest. The estimated error in this quantity, E Q,est , and the exact error E Q are plotted for combinations respectively. As expected it is not until the basis include five SD modes that the benefit from the higher eigenvalue gives an impact on the estimate. The result from the adaptive procedure is plotted in Fig. 22. As the figure shows, this example is challenging for the adaptive strategy, which, in fact, performs worse than the pure POD basis. This can be explained by the eigenvalue spectrum (Fig. 2) which has multiple eigenvalues with the same magnitude. Even though some SD modes are added at steps 7, 10, 14, and 15, when effectiveness of the POD basis decreases, they do not decrease the estimate much since the eigenvalue stays the same.

Conclusions and outlook
In this paper we have presented an a posteriori error estimate for the error introduced by NMR in the context of homog-enization of porous media. In particular, we proposed and investigated a combined basis, consisting of SD modes and POD modes, as an extension to the reduced model presented in Jänicke et al. [12]. The explicit error estimator proposed in Ekre et al. [24] is presented for the combined basis for estimation of the NMR error. The combined basis utilizes both the "residual minimizing" property of the POD modes and the "estimator improving" property of the SD modes. It is the estimated error that will be considered in the pertinent error control, i.e., it is the estimated error that will be compared to the pre-set tolerance. Therefore, the best procedure is the one that shows the most efficient convergence in the (guaranteed) error estimator. A procedure for adaptively choosing the "correct" next mode is presented. This procedure indicates whether the next mode should come from the pool of SD modes or POD modes. The performance of the estimator, with both (i) energy norm and (ii) user-defined quantities of interest as the target, was demonstrated by examples. We show that it is possi-  ble to find an improved composition of modes in order to minimize the estimated error. The results also show that the proposed adaptive mode selection strategy works well, but has problems when there are multiple eigenvalues close in the spectrum. For future work it would be of interest to utilize the error estimate, and the adaptive scheme, for the full two-scale FE 2 problem.
Funding Open access funding provided by Chalmers University of Technology.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A NMR error estimate for generic error equation
Consider the following generic error equation; Find χ ∈P s.t.Â where the generic residual is defined as for generic functions M t (x, t) ∈ P and M 0 (x), M T (x) ∈ P . The identity R (q) = R ( C q) follows from (i) the assumption 8 that R (q R ) = 0 ∀q R ∈ P ,R , and (ii) the definition of the projection operator. In order to estimate χ we need to estimate χ a , χ(•, 0) m , and χ(•, 0) m in terms of the data in the right hand side, i.e. in terms of M t , M 0 , and M T . We note that sinceÂ(•, •) localizes in time we can express the generic error equation explicitly for t = 0, t ∈ I and t = T and obtain explicit computable bounds for the error.

A.1 Generic error equation at t = 0 and t = T
At the endpoints of the time interval the generic error equation reduced to the following two equations, with τ = 0 and τ = T . Find χ(τ )=:χ τ ∈ P s.t. = 2m ( C M τ , χ τ ) from which we obtain an upper estimate for the contributions to χ , viz.

A.2 Generic error equation at t ∈ I
Inside the time interval the generic error equation reduced to the following. Find χ(t)=:χ t ∈ P s.t.
from which we obtain the following upper bound using Cauchy-Schwarz inequality In order to relate • a to • m we utilize the properties of the spectral basis, and write χ t = N a=1 ξ a ϕ a where ξ a = ξ a (t) are coefficients, and where ϕ a = ϕ a (x) are the eigenmodes from the solution to the following eigenvalue problem where λ a are the eiganvalues in increasing order (λ 1 ≤ λ 2 . . .) and ϕ a the corresponding eigenmodes. From the generic error Eq. (78) and the Galerkin-like property of the residual (cf. Sects. 4.2 and 4.3), we note that ξ a ≡ 0, a = 1, 2, . . . , N SD R , and consequently we note that From the relationship between a (•, •) and m (•, •) given by the eigenvalue problem in (80) and the properties of the projection we obtain the following upper estimate: Remark We note that the estimate in (82) can be improved by replacing λ N SD R with λ N SD R +1 , i.e. with a larger eigenvalue. However, in practice we will only consider the N SD R eigenvalues that are already computed as part of constructing the reduced basis.
Equations (79) and (82) now results in the following estimate of χ t a χ t a ≤ C M t m λ N SD R . (83)

A.3 Final NMR error estimate of the generic error
We maybe now combine the definition of the norm in Eq. (52) with the partial estimates from Eqs. (77) and (83)  (84)