Localized bases for finite dimensional homogenization approximations with non-separated scales and high-contrast.

Localized bases for finite dimensional homogenization approximations with non-separated scales and high-contrast.

Abstract

We construct finite-dimensional approximations of solution spaces of divergence form operators with -coefficients. Our method does not rely on concepts of ergodicity or scale-separation, but on the property that the solution space of these operators is compactly embedded in if source terms are in the unit ball of instead of the unit ball of . Approximation spaces are generated by solving elliptic PDEs on localized sub-domains with source terms corresponding to approximation bases for . The -error estimates show that -dimensional spaces with basis elements localized to sub-domains of diameter (with ) result in an accuracy for elliptic, parabolic and hyperbolic problems. For high-contrast media, the accuracy of the method is preserved provided that localized sub-domains contain buffer zones of width where the contrast of the medium remains bounded. The proposed method can naturally be generalized to vectorial equations (such as elasto-dynamics).

1Introduction

Consider the partial differential equation

where is a bounded subset of with a smooth boundary (e.g., ) and is symmetric and uniformly elliptic on . It follows that the eigenvalues of are uniformly bounded from below and above by two strictly positive constants, denoted by and . Precisely, for all and ,

In this paper, we are interested in the homogenization of (and its parabolic and hyperbolic analogues in Sections Section 4 and Section 5), but not in the classical sense, i.e., that of asymptotic analysis [9] or that of or -convergence ([47], [57]) in which one considers a sequence of operators and seeks to characterize limits of solution. We are interested in the homogenization of in the sense of “numerical homogenization,” i.e., that of the approximation of the solution space of by a finite-dimensional space.

This approximation is not based on concepts of scale separation and/or of ergodicity but on compactness properties, i.e., the fact that the unit ball of the solution space is compactly embedded into if source terms () are integrable enough. This higher integrability condition on is necessary because if spans , then the solution space of is (and it is not possible to obtain a finite dimensional approximation subspace of with arbitrary accuracy in -norm). However, if spans the unit ball of , then the solution space of shrinks to a compact subset of that can be approximated to an arbitrary accuracy in -norm by finite-dimensional spaces [10] (observe that if , then the solution space is a closed bounded subset of , which is known to be compactly embedded into ).

The identification of localized bases spanning accurate approximation spaces relies on a transfer property obtained in [10]. For the sake of completeness, we will give a short reminder of that property in Section 2. In Section 3, we will construct localized approximation bases with rigorous error estimates (under no further assumptions on than those given above). In sub-Section 3.4, we will also address the high-contrast scenario in which is allowed to be large. In Sections Section 4 and Section 5, we will show that the approximation spaces obtained by solving localized elliptic PDEs remain accurate for parabolic and hyperbolic time-dependent problems. We refer to Section 6 for numerical experiments. We refer to Section Appendix B of the Appendix for further discussion and a proof of the strong compactness of the solution space when the range of is a closed bounded subset of with (this notion of strong compactness constitutes a simple but fundamental link between classical homogenization, numerical homogenization and reduced order modeling).

2A reminder on the flux-norm and the transfer property.

Recall that the key element in and convergence is a notion of “compactness by compensation” combined with convergence of fluxes. Here, the notion of compactness is combined with a flux-norm introduced in [10].

The flux-norm. We will now give a short reminder on the flux-norm and its properties.

The following proposition shows that the flux-norm is equivalent to the energy norm if and .

Motivations behind the flux-norm: There are three main motivations behind the introduction of the flux norm.

• The flux-norm allows to obtain approximation error estimates independent from both the minimum and maximum eigenvalues of . In fact, the flux-norm of the solution of is independent from altogether since

• The in the -norm is explained by the fact that in practice, we are interested in fluxes (of heat, stress, oil, pollutant) entering or exiting a given domain. Furthermore, for a vector field , , which means that the flux entering or exiting is determined by the potential part of the vector field.

• Classical homogenization is associated with two types of convergence: convergence of energies (-convergence [33]) and convergence of fluxes ( or -convergence [47]). Similarly, one can define an energy norm and a flux-norm.

The transfer property. For , a finite dimensional linear subspace of , we define

Note that is a finite dimensional subspace of .

The usefulness of can be illustrated by considering so that . Then, and therefore can be chosen as, e.g., the standard piecewise linear FEM space, on a regular triangulation of of resolution , with nodal basis . The space is then defined by its basis determined by

Equation shows that the approximation error estimate associated with the space and the problem with arbitrarily rough coefficients is (in -flux norm) equal to the approximation error estimate associated with piecewise linear elements and the space . More precisely,

where does not depend on .

We refer to [22], [25] and [11] for recent results on finite element methods for high contrast () but non-degenerate () media under specific assumptions on the morphology of the (high-contrast) inclusions (in [22], the mesh has to be adapted to the morphology of the inclusions). Observe that the proposed method remains accurate if the medium is both of high contrast and degenerate (), without any further limitations on , at the cost of solving PDEs over the whole domain .

3Localization of the transfer property.

The elliptic PDEs have to be solved on the whole domain . Is it possible to localize the computation of the basis elements to a neighborhood of the support of the elements ? Observe that the support of each is contained in a ball of center (the node of the coarse mesh associated with ) of radius . Let . Solving the PDEs on sub-domains of (containing the support of ) may, a priori, increase the error estimate in the right hand side of . This increase can, in fact, be linked to the decay of the Green’s function of the operator . The slower the decay, the larger the degradation of those approximation error estimates. Inspired by the strategy used in [35] for controlling cell resonance errors in the computation of the effective conductivity of periodic or stochastic homogenization (see also [36]), we will replace the operator by the operator in the left hand side of in order to artificially introduce an exponential decay in the Green’s function. A fine tuning of is required because although a decrease in improves the decay of the Green function, it also deteriorates the accuracy of the transfer property. In order to limit this deterioration, we will transfer a vector space with a higher approximation order than the one associated with piecewise linear elements. Let us now give the main result.

3.1Localized bases functions.

Let . Let be an approximation sub-vector space of such that

• is spanned by basis functions (with ) with supports in where, the are the nodes of a regular triangulation of of resolution .

• satisfies the following approximation properties: For all

and for all

• For all ,

• For all coefficients ,

Through this paper, we will write any constant that does not depend on (but may depend on , , and the essential supremum and infimum of the maximum and minimum eigenvalues of over ). Let and . For each basis element of let be the solution of

Let

be the linear space spanned by the elements .

On Localization. We refer to [22], [25] and [6] for recent localization results for divergence-form elliptic PDEs. The strategy of [22] is to construct triangulations and finite element bases that are adapted to the shape of high conductivity inclusions via coefficient dependent boundary conditions for the subgrid problems (assuming to be piecewise constant and the number of inclusions bounded). The strategy of [25] is to solve local eigenvalue problems, observing that only a few eigenvectors are sufficient to obtain a good pre-conditioner. Both [22] and [25] require specific assumptions on the morphology and number of inclusions. The idea of the strategy is to observe that if is piecewise constant and the number of inclusions bounded, then is locally away from the interfaces of the inclusions. The inclusions can then be taken care of by adapting the mesh and the boundary values of localized problems or by observing that those inclusions will affect only a finite number of eigenvectors.

The strategy of [6] is to construct Generalized Finite Elements by partitioning the computational domain into to a collection of preselected subsets and compute optimal local bases (using the concept of -widths [55]) for the approximation of harmonic functions. Local bases are constructed by solving local eigenvalue problems (corresponding to computing eigenvectors of where is the restriction of -harmonic functions from onto , is the adjoint of , and is a sub-domain of surrounded by a larger sub-domain ). The method proposed in [6] achieves a near exponential convergence rate (in the number of pre-computed bases functions) for harmonic functions. Non-zero right hand sides () are then taken care of by solving (for each different ) particular solutions on preselected subsets with a constant Neumann boundary condition (determined according to the consistency condition).

As explained in Remark ?, the near exponential rate of convergence observed in [6] is explained by the fact that the source space considered in [6] is more regular than (since [6] requires the computation particular (local) solutions for each right hand sides and each non-zero boundary conditions, the basis obtained in [6] is in fact adapted to -harmonic functions away from the boundary). The strategy proposed here can also be used to achieve exponential convergence for analytic source terms by employing higher-order basis functions in . Furthermore, as shown in sections Section 4, Section 5 and Section 3.4 the method proposed here allows for the numerical homogenization of time-dependent problems (because it does not require the computation of particular solutions for different source or boundary terms) and can be extended to high-contrast media. We also note that the basis functions are simpler and cheaper to compute (equation ) than the eigenvectors of required by [6]. We refer to page 16 of [6] for a discussion on the cost of this added complexity.

3.2On Numerical Homogenization.

By now, the field of numerical homogenization has become large enough that it is not possible to give an exhaustive review in this short paper. Therefore, we will restrict our attention to works directly related to our work.

- The multi-scale finite element method [40] can be seen as a numerical generalization of this idea of oscillating test functions found in -convergence. A convergence analysis for periodic media revealed a resonance error introduced by the microscopic boundary condition [40]. An over-sampling technique was proposed to reduce the resonance error [40].

- Harmonic coordinates play an important role in various homogenization approaches, both theoretical and numerical. These coordinates were introduced in [42] in the context of random homogenization. Next, harmonic coordinates have been used in one-dimensional and quasi-one-dimensional divergence form elliptic problems [?], allowing for efficient finite dimensional approximations. The connection of these coordinates with classical homogenization is made explicit in [2] in the context of multi-scale finite element methods. The idea of using particular solutions in numerical homogenization to approximate the solution space of appears to have been first proposed in reservoir modeling in the 1980s [16], [61] (in which a global scale-up method was introduced based on generic flow solutions i.e., flows calculated from generic boundary conditions). Its rigorous mathematical analysis was done only recently [49] and is based on the fact that solutions are in fact -regular with respect to harmonic coordinates (recall that they are -regular with respect to Euclidean coordinates). The main message here is that if the right hand side of is in , then solutions can be approximated at small scales (in -norm) by linear combinations of (linearly independent) particular solutions ( being the dimension of the space). In that sense, harmonic coordinates are only good candidates for being linearly independent particular solutions.

The idea of a global change of coordinates analogous to harmonic coordinates has been implemented numerically in order to up-scale porous media flows [27]. We refer, in particular, to a recent review article [16] for an overview of some main challenges in reservoir modeling and a description of global scale-up strategies based on generic flows.

- In [24], the structure of the medium is numerically decomposed into a micro-scale and a macro-scale (meso-scale) and solutions of cell problems are computed on the micro-scale, providing local homogenized matrices that are transferred (up-scaled) to the macro-scale grid. This procedure allows one to obtain rigorous homogenization results with controlled error estimates for non-periodic media of the form (where is assumed to be smooth in and periodic or ergodic with specific mixing properties in ). Moreover, it is shown that the numerical algorithms associated with HMM and MsFEM can be implemented for a class of coefficients that is much broader than . We refer to [34] for convergence results on the Heterogeneous Multiscale Method in the framework of and -convergence.

- More recent work includes an adaptive projection based method [48], which is consistent with homogenization when there is scale separation, leading to adaptive algorithms for solving problems with no clear scale separation; fast and sparse chaos approximations of elliptic problems with stochastic coefficients [60]; finite difference approximations of fully nonlinear, uniformly elliptic PDEs with Lipschitz continuous viscosity solutions [19] and operator splitting methods [4].

- We refer to [13] (and references therein) for most recent results on homogenization of scalar divergence-form elliptic operators with stochastic coefficients. Here, the stochastic coefficients are obtained from stochastic deformations (using random diffeomorphisms) of the periodic and stationary ergodic setting.

3.3Proof of Theorem .

For each basis element of , let be the solution of

The following Proposition will allow us to control the impact of the introduction of the term in the transfer property. Observe that the domain of PDE is still (our next step will be to localize it to ).

Let . We have

Define to be the energy norm . Multiplying by and integrating by parts, we obtain that

Write and let and be the solutions of and with Dirichlet boundary conditions on . Then, we obtain by integration by parts and the Cauchy-Schwartz inequality that

Using , we can choose so that

Using , we can choose so that

we conclude the proof of the approximation by observing that . Let us now prove Equation . First, observe that Equation and the triangular inequality imply that

Next, we obtain from and Poincaré inequality and

and

We conclude by combining equations and with .

We will now control the error induced by the localization of the elliptic problem . To this end, for each each basis element of write the intersection of the support of with and let be a subset of containing such that . Let also be the solution of

For , write the Euclidean distance between the sets and .

We refer to Section Appendix A of the Appendix for the proof of Proposition ?.

Taking (we use the particular notation because our proof of accuracy requires that specific constant to be large enough, i.e., larger than a constant depending on the parameter appearing in the right hand side of and the parameter describing the balls containing the support of the basis functions introduced in subSection 3.1) and in equation of Proposition ?, we obtain for large enough (but independent from ) that

Let be the solution of in . Using Proposition ?, we obtain that there exist coefficients such that

and

Using the triangle inequality, it follows that

whence, from Cauchy-Schwartz inequality,

Combining with , we obtain that

Using in , we obtain that

Observe that it is the exponential decay in that allows us to compensate for the large term on the right hand side of via . This concludes the proof of Theorem ?.

3.4On localization with high-contrast.

The constant in the approximation error estimate depends, a priori, on the contrast of . Is it possible to localize the computation of bases for when the contrast of is high? The purpose of this subsection is to show that the answer is yes provided that there is a buffer zone between the boundaries of localization sub-domains and the supports of the elements where the contrast of remains bounded. More precisely, assume that is the disjoint union of and . Assume that holds only on , and that on we have

where can be arbitrarily large. Practical examples include media characterized by a bounded contrast background with high conductivity inclusions or channels. Let be the solution of

Let

be the linear space spanned by the elements . For each , define to be the largest number such that there exists a subset such that: the closure of contains the support of , is a subset of (where are the set of points of that are at distance at most for ), and is a subset of . If no such subset exists we set . can be interpreted as the non-high-contrast buffer distance between the support of and the boundary of . We refer to Figure ? for illustrations of the buffer distance.

The proof of Theorem ? is similar to that of Theorem ?, but it requires a precise tracking of the constants involved. Because of the close similarity we will not include the proof in this paper but only give its main lines. First, the proof of Proposition ? remains unchanged as the constants in and do not depend on the maximum eigenvalue of the conductivity . Only the proof of Proposition ? has to be adapted and the part of the proof below Proposition ? remains unchanged. This requires an application of the elements of lemmas ?, ?, ? and ? to buffer sub-domains . The main point is to observe that the decay of the Green’s function in can be bounded independently from (due to the maximum principle).

Observe that the choice of the sub-domain in can be chosen to be the same as in if its intersection with high contrast inclusions is void (i.e., if the maximum eigenvalue of over remains bounded independently from ); otherwise the choice of in has to be enlarged (when compared to that associated with ) to contain the high-contrast inclusion (plus its buffer).

4The basis remains accurate for parabolic PDEs.

The computational gain of the method proposed in this paper is particularly significant for time-dependent problems. One such problem is the parabolic equation associated with the operator . More precisely, consider the time-dependent partial differential equation

where and satisfy the same assumptions as those associated with PDE , for some and .

Let be the finite-dimensional approximation space defined in . Let be the finite element solution of , i.e., can be decomposed as

and solves for all

with . Write

The proof is a generalization of the proof found in [50] (in which approximation spaces are constructed via harmonic coordinates). Let be the bilinear form on defined by

Observe that for all ,

Writing , we deduce that for ,

Using in and integrating, we obtain that

Using Minkowski’s inequality, we deduce that

Similarly,

Using Cauchy-Schwartz and Minkowski inequalities in , we obtain that

Take to be the projection of onto with respect to the bilinear form . Observing that with , we obtain from Theorem ? that

Let us now show (using a standard duality argument) that

Choose to be the solution of the following linear problem: For all

Taking in , we obtain that

Hence by Cauchy Schwartz inequality and ,

Using Theorem ? again, we obtain that

which concludes the proof of Theorem ?.

Discretization in time. Let be a discretization of with time-steps . Write , the subspace of , such that

Write , the solution in of the following system of implicit weak formulation (such that ): For each and ,

Then, we have the following theorem

The proof of Theorem ? is similar to that of Theorem 1.6 of [50] and will not be given here. Observe that homogenization in space allows for a discretization in time with time steps without compromising the accuracy of the method.

5The basis remains accurate for hyperbolic PDEs.

Consider the hyperbolic partial differential equation

where , , and are defined as in Section 4. In particular, is assumed to be only uniformly elliptic and bounded (). We will further assume that is uniformly bounded from below and above ( and ). It is straightforward to extend the results presented here to nonzero boundary conditions (provided that frequencies larger than remain weakly excited, because the waves equation preserves energy and homogenization schemes can not recover energies put into high frequencies, see [51]). For the sake of conciseness, we will give those results with zero boundary conditions. PDE corresponds to acoustic wave equations in a medium with density and bulk modulus .

Let be the finite-dimensional approximation space defined in . Let be the finite element solution of , i.e., can be decomposed as

and solves for all

where

Let be the bilinear form on defined in . Observe that for all ,

Taking as a test function in and integrating in time, we deduce that for ,

where . Taking the derivative of the hyperbolic equation for in time, we obtain that

Integrating against the test function and observing that , we also obtain that

Similarly, we obtain that

Take to be the projection of onto with respect to the bilinear form . Observing that with , we obtain from and Theorem ? that

Furthermore, using the same duality argument as in the parabolic case, we obtain that

Using Cauchy-Schwartz and Minkowski inequalities and the above estimates in , we obtain that

We conclude using Grownwall’s lemma.

6Numerical experiments.

6.1Elliptic equation.

We compute the solutions of up to time on the fine mesh and in the finite-dimensional approximation space defined in . The physical domain is the square . Global equations are solved on a fine triangulation with nodes and triangles.

The elements of sub-Section 3.1 are weighted extended B-splines (WEB) [38] (obtained by tensorizing one-dimensional elements and using weight function to enforce the Dirichlet boundary condition). The order of accuracy is not affected by the choice of weight function given that the boundary is piecewise smooth. Our motivation for using WEB elements lies in the fact that, with those elements, (Dirichlet) boundary conditions become simple to enforce. This being said, any finite elements satisfying the properties , , and would be adequate [14].

We write the size of the coarse mesh. Elements are obtained by solving on localized sub-domains of size . Table ? shows errors with for given by (Example 1 of Section 3 of [49], trigonometric multi-scale, see also [45]), i.e., for