An offline-online strategy for multiscale problems with random defects
AAn offline-online strategy for multiscaleproblems with random defects ∗ Axel M˚alqvist † Barbara Verf¨urth ‡ Abstract.
In this paper, we propose an offline-online strategy based on the Localized Orthogonal Decompo-sition (LOD) method for elliptic multiscale problems with randomly perturbed diffusion coefficient. We consider aperiodic deterministic coefficient with local defects that occur with probability p . The offline phase pre-computesentries to global LOD stiffness matrices on a single reference element (exploiting the periodicity) for a selection ofdefect configurations. Given a sample of the perturbed diffusion the corresponding LOD stiffness matrix is thencomputed by taking linear combinations of the pre-computed entries, in the online phase. Our computable errorestimates show that this yields a good coarse-scale approximation of the solution for small p . Moreover, extensivenumerical experiments illustrate that relative errors of a few percent are achieved up to at least p = 0 .
1. Thismakes the proposed technique attractive already for moderate sample sizes in a Monte Carlo simulation.
Key words. numerical homogenization, multiscale method, finite elements, random perturbations
AMS subject classifications.
1. Introduction
Many modern materials include some fine composite structure to achieve enhanced properties.Examples include fiber reinforced structures in mechanics as well as mechanical, acoustic or op-tical metamaterials. The materials are often highly structured, but mistakes in the fabricationprocess lead to defects. A major question is the robustness of the desired material properties un-der such defects. Mathematically speaking, we are interested in the solution of partial differentialequations (PDEs) with multiscale, randomly perturbed coefficients.In this paper, we study the following elliptic multiscale problem: Find u : D → R such that − ∇ · ( A ( x ) ∇ u ( x )) = f ( x ) in D (1.1)with suitable boundary conditions. Here, D is a spatial domain in R d and f ∈ L ( D ). Themultiscale coefficient A ∈ L ∞ ( D, R ) is a sample of a randomly perturbed coefficient. Moreprecisely we assume that A is a realization of the form A ( x, ω ) = A ε ( x ) + b p ( x, ω ) B ε ( x ) , (1.2) ∗ This work was initiated while BV enjoyed the kind hospitality of Chalmers University, Gothenburg. AM isfunded by the Swedish Research Council and the G¨oran Gustafsson foundation for Research in Natural Sciencesand Medicine. BV is funded by the German Research Foundation (DFG) – Project-ID 258734477 – SFB 1173and Klaus-Tschira foundation as well as by the Federal Ministry of Education and Research (BMBF) and theBaden-W¨urttemberg Ministry of Science as part of the Excellence Strategy of the German Federal and StateGovernments. † Department of Mathematical Sciences, Chalmers University of Technology and University of Gothenburg, 41296 G¨oteborg, Sweden ‡ Institut f¨ur Angewandte und Numerische Mathematik, Karlsruher Institut f¨ur Technologie, Englerstr. 2, 76131Karlsruhe, Germany a r X i v : . [ m a t h . NA ] F e b igure 1.1.: Two examples for weakly random coefficients: random checkerboard ( ε = 2 − , p =0 .
1, left) and periodic inclusions with random “erasure” ( ε = 2 − , p = 0 .
1, right)where A ε , B ε are deterministic multiscale coefficients and b p ( x, · ) is a Bernoulli law with prob-ability p , cf. [3]. Detailed assumptions on the problem data and the form of A are given inSection 2 below. Two important examples of this setup are illustrated in Figure 1.1. On the left, A is generated from a constant A ε by introducing square spots with length ε and probability p .On the right, A ε is made of a background value and periodic square inclusion repeating witha periodicity length ε . In this case, A is generated by randomly setting some of the inclusionsto the background value, thus “erasing” them. The considered model of so-called weakly ran-dom coefficients, characterized by small values of p , also covers other defect possibilities of theinclusions like a change of value, a (fixed) shift or a (fixed) change of the geometry.In the context of materials with defects, one is interested in extracting statistical informationabout the solution u . There are many different uncertainty quantification techniques for PDEswith random coefficients, see, e.g., [4, 15, 24] for overviews. In the following, we focus on MonteCarlo (MC)-type approaches such as Quasi Monte Carlo or Multilevel Monte Carlo (MLMC)[5, 7, 29, 10]. Hence, we are interested in (approximate) solutions to (1.1) for many samples (i.e.,realizations) of A . Due to the multiscale nature of A , standard discretization schemes like thefinite element method would require the mesh to resolve all fine-scale features. Consequently,already the computation of a few solutions to (1.1) becomes prohibitively costly. In contrast,computational multiscale methods such as the Localized Orthogonal Decomposition (LOD) [20,25, 26, 1] yield faithful coarse-scale approximations with feasible effort after pre-computationof a generalized finite element basis. However, these basis functions incorporate knowledgeabout the multiscale coefficient A and thus, need be constructed anew for each realization ingeneral. This makes it difficult to apply multiscale methods to stochastic problems. Recentlythere have been several attempts to circumvent this difficulty, for example the combination of theMultiscale Finite Element Method with Multilevel Monte Carlo [8] or low-rank approximation[28], the multiscale data-driven stochastic method [30], an approach for a quasi-local homogenizedcoefficient [13, 14], and a sparse compression of the expected solution operator [12]. In the contextof the aforementioned LOD, the recent works [18, 17] deal with rare defects. They propose tocompute the multiscale basis for the unperturbed deterministic coefficient and to update thisbasis only locally for each particular sample. More precisely, given a sample coefficient A , acomputable error indicator shows whether the pre-computed basis is sufficiently good or if a new(improved) basis needs to be computed for that particular sample. If p is small enough this2echnique becomes competitive. A pre-computed, deterministic multiscale basis is also the keyidea of the Multiscale Finite Element approach of [22]. An asymptotic expansion in the randomvariable is used for the numerical analysis of [22] as well as to approximate the effective coefficientof stochastic homogenization [2, 3, 23, 21] or to reduce its variance [6, 21].This contribution is an attempt to make multiscale methods useful for a wider range of randomproblems. Instead of having one reference coefficient, as in [17], we build a “basis” of referencecoefficients { A i } Ni =0 and pre-compute and store the corresponding LOD basis functions { λ −C ( A i ) } Ni =0 , λ being the finite element basis function and C ( B ) the LOD correction based oncoefficient B . This is done in the offline phase. We consider problems with periodic structureso that the same set of precomputed basis functions can be used in the entire domain. Givensamples of the form A = (cid:80) Ni =0 µ i A i we let (cid:80) Ni =0 µ i C ( A i ) λ approximate C ( A ) λ in the onlinephase. This allows for rapid assembly of the LOD stiffness matrices and thereby solution of theproblem given samples A . We demonstrate theoretically that the error C ( A ) λ − (cid:80) Ni =0 µ i C ( A i ) λ is small for small p in one dimension and provide a computable error indicator of this quantity inhigher dimensions. We also present numerical experiments for a large variety of configurationsof the diffusions which shows relative root mean square errors up to 3% for defect probabilitiesof p = 0 . A . In Section 3, we review the Petrov-Galerkin Localized Orthogonal Decomposition(PG-LOD) and introduce our new offline-online strategy. A priori error estimates for the newmethod are presented in Section 4. We discuss several implementation details with a focus oncomputational efficiency in Section 5. Extensive numerical experiments in Section 6 showcasethe attractive properties of the method and also illustrate our theoretical findings.
2. Problem formulation
In the following, we detail the setting associated with (1.1). We first pose the problem fora fixed event ω in a probability space Ω and then discuss the specific form of randomness inthe coefficient. By slight abuse of notation, we will omit the random variable in the followingexposition for a fixed, but arbitrary sample. For simplicity we let D = [0 , d ⊂ R d be the unit cell. We assume that f ∈ L ( D ) and that therealization A ∈ L ∞ ( D, R ) is uniformly bounded and elliptic, i.e.,0 < α := ess inf x ∈ D A ( x ) , ∞ > β := ess sup x ∈ D A ( x ) . (2.1)We introduce a function space V where we seek a solution of (1.1) in weak form. In this paperwe mainly consider a conforming finite element space V := V h ⊂ H , ( D ) = { v ∈ H ( D ) | v is periodic , ˆ D v = 0 } defined on a computational mesh T h which resolves the variations in A . Further, we assume T h to be a conforming (i.e., without hanging nodes and edges) and shape regular quadrilateralmesh that can additionally be wrapped into a mesh on the torus without hanging nodes oredges. However, it is also possible to choose V = H , ( D ) and the following analysis will still3o through. The weak form of (1.1) reads as follows: find u ∈ V such that a ( u, v ) = F ( v ) for all v ∈ V, (2.2)where a ( u, v ) := ˆ D A ( x ) ∇ u ( x ) · ∇ v ( x ) dx, F ( v ) := ˆ D f ( x ) v ( x ) dx. Due to the constraint ´ D u ( x, ω ) dx = 0, existence and uniqueness of a solution u is guaranteedby the Lax-Milgram lemma. We will frequently use the energy norm (cid:107) · (cid:107) A := a ( · , · ) / and itsrestriction (cid:107) · (cid:107) A,S to a subdomain S ⊂ D in the following. Further, ( · , · ) S denotes the usual L -scalar product on S where we omit the subscript if S = D . Remark 2.1.
We consider periodic boundary conditions and box-type domains in this paperto fully exploit the underlying structure in A . However, with additional computational effort,Dirichlet or Neumann boundary conditions as well as more general Lipschitz domains can betreated, see Remark 3.1. In particular, the error analysis is not restricted to periodic boundaryconditions or box-type domains. As mentioned above, we are interested in solving (2.2) for many different (random) choices of A .We now give more details on the form (1.2) of A ( x, ω ), similar to the “weakly” random settingof [3]. We assume that A ε and B ε are deterministic multiscale coefficients. More specifically, welet A ε ( x ) = A per ( x/ε ) where A per is 1-periodic and we assume that ε = 1 /n with n ∈ N , n (cid:29) B ε ( x ) = B per ( x/ε ) is assumed for B ε . The periodicity assumption – togetherwith the box-type domain – is imposed for efficiency reasons of our method, where we emphasizethat generalizations are possible, see Remark 3.2.The deterministic coefficients are assumed to satisfy spectral bounds similar as (2.1), i.e.,0 < α ≤ ess inf x ∈ D A per ( x ) , ∞ > β ≥ ess sup x ∈ D A per ( x ) . (2.3)and 0 < α ≤ ess inf x ∈ D ( A per ( x ) + B per ( x )) , ∞ > β ≥ ess sup x ∈ D ( A per ( x ) + B per ( x )) . (2.4)The random character of A ( x, ω ) in (1.2) is encoded in b p ( x, ω ) for which we assume b p ( x, ω ) = (cid:88) j ∈ I χ ε ( j + Q ) ( x )ˆ b jp ( ω ) . (2.5)Here, χ denotes the characteristic function, Q ⊆ [0 , d and I := { k ∈ Z d | ε ( k + Q ) ⊂ D } . Finally,ˆ b jp are independent random variables adhering to a Bernoulli distribution with probability p ,i.e., ˆ b jp = 0 with probability 1 − p and ˆ b jp = 1 with probability p . Clearly, for p →
0, thedefects/perturbations encoded in b p become rare events. This choice of b p together with theassumptions (2.3)–(2.4) guarantee that each realization A ( x, ω ) satisfies (2.1). We close thissection by giving two examples for this setting, namely the formal definition of the coefficientsdepicted in Figure 1.1 above. Example 2.2 (Random checkerboard) . Recall that the coefficient in Figure 1.1, left, is piece-wise constant on a square mesh T ε . On each square element, the value of A is picked randomly aseither α with probability 1 − p or as β with probability p . This can be described in the form (1.2)and (2.5) with Q = [0 , d , A per = α and B per = β − α .4 xample 2.3 (Periodic coefficient with random defects) . Realizations as depicted in Figure 1.1,right, can be formalized in the following way. We define A per : [0 , d → R via A per ( y ) := (cid:40) β y ∈ [0 . , . d ,α else . Further, we pick Q = [0 . , . d and B per = α − β . Note that B per is only added in theshifted and scaled copies of Q . Clearly, any other value ˜ β ∈ [ α, β ] as defect can be modeled bythe choice B per = ˜ β − α . Even a value 0 < ˜ β / ∈ [ α, β ] is possible for the defects by changingthe spectral bounds. If the defect changes the shape of the inclusion, we have to define Q and B per accordingly. For example, imagine that a defect means that the value β is taken in (scaledand shifted copies of) [0 . , d . Then, A per is left unchanged, we set Q = [0 , d and define B per : Y → R via B per ( y ) := α − β y ∈ [0 . , . d ,β − α y ∈ [0 . , d , .
3. Offline-online strategy for the PG-LOD
In this section, we first review the Petrov-Galerkin Localized Orthogonal Decomposition (PG-LOD) in Sections 3.1 and 3.2 and then present our new offline-online strategy in Section 3.3.Throughout this paper, we further use the notation a (cid:46) b if a ≤ cb with a generic constant c that only depends on the shape regularity of the mesh, the domain D , or the space dimension d . Let T H be a coarse, shape regular, quasi-uniform and conforming quadrilateral mesh of thedomain D . We further assume that T H can be wrapped into a conforming mesh of the torus,i.e., no hanging nodes and edges occur over the periodic boundary. Let H = max T ∈T H diam T denote the mesh size. The standard lowest-order finite element space on T H is given as V H := H , ( D ) ∩ Q ( T H ) , where Q ( T H ) denotes the space of T H -piecewise polynomials of coordinate degree at most 1.We further assume that the fine mesh T h is a refinement of T H such that the finite element spacesare nested as V H ⊂ V h .We further introduce a notion of element patches. For an arbitrary subdomain S ⊂ D and m ∈ N , we define patches U m ( S ) ⊂ D inductively as U ( S ) = S, U m +1 ( D ) = (cid:91) { T ∈ T H | U m ( S ) ∩ T (cid:54) = ∅} . In this definition, T H is interpreted as a mesh of the torus such that patches are continued overthe periodic boundary [27]. For S = T with T ∈ T H we call U m ( T ) the m -layer element patchand we refer to Figure 3.1 for a visualization. By the quasi-uniformity of T H we further notethat max T ∈T H card { K ∈ T H | K ∈ U m ( T ) } (cid:46) m d . (3.1)Fine-scale features are not captured in the coarse space V H , i.e., a standard FEM on the coarsescale applied to (2.2) does not yield faithful approximations. We will characterize fine-scale parts5igure 3.1.: A mesh element T = U ( T ) (dark blue) and its patches U ( T ) (intermediate blue)and U ( T ) (light blue). The light blue squares on the very right belong to U ( T )because of the periodic continuation.of functions in V as the kernel of a suitable interpolation operator. We now introduce the requiredproperties as well as an appropriate example. Let I H : V → V H denote a bounded local linearprojection operator, i.e., I H ◦ I H = I H , with the following stability and approximation propertiesfor all v ∈ V H − (cid:107) v − I H v (cid:107) L ( T ) + (cid:107)∇ I H v (cid:107) L ( T ) (cid:46) (cid:107)∇ v (cid:107) L ( U ( T )) . (3.2)A possible choice (which we use in our implementation of the method) is to define I H := E H ◦ Π H ,where Π H : V → Q ( T H ) denotes the (local) L -projection. E H is the averaging operator thatmaps discontinuous functions in Q ( T H ) to V H by assigning to each free vertex the arithmeticmean of the corresponding function values of the neighboring cells, that is, for any v ∈ Q ( T H )and any vertex z of T H ,( E H ( v ))( z ) = (cid:88) T ∈T H , z ∈ T v | T ( z ) (cid:30) card { K ∈ T H , z ∈ K } . Note that again T H is understood as a mesh of the torus. For further details on suitable inter-polation operators we refer to [11]. Denote by V f = ker I H the kernel of I H which characterizes the functions with possible fine-scalevariations in V . Further, we introduce the following straightforward restriction of V f to patches U m ( S ) via V f ( U m ( S )) := { v ∈ V f | v | D \ U m ( S ) = 0 } . The LOD is based on (truncated) correction operators C m ( A ) defined by C m ( A ) v = (cid:88) T ∈T H C m,T ( A ) v, where the so-called element correction operator C m,T ( A ) : V → V f ( U m ( T )) associated with thecoefficient A solves (cid:0) A ∇ ( C m,T ( A ) v ) , ∇ v f (cid:1) U m ( T ) = (cid:0) A ∇ v, ∇ v f (cid:1) T for all v f ∈ V f ( U m ( S )) . (3.3)6ote that these local problems are well-posed by the Lax-Milgram lemma due to the uniformellipticity of A . The multiscale space V ms H,m is now constructed as V ms H,m := V H − C m ( A ) V H . Denoting by N the set of vertices of T H (understood as a mesh of the torus) and { λ z } z ∈N thenodal basis of V H , { λ z − C m ( A ) λ z } z ∈N is a basis of V ms H,m .The Petrov-Galerkin LOD (PG-LOD) for (2.2) now reads as: Find u ms m ∈ V ms H,m such that a ( u ms m , v ) = F ( v ) for all v ∈ V H , (3.4)where we can write u ms m = u Hm − C m ( A ) u Hm with u Hm ∈ V H . In this Petrov-Galerkin variant onlythe ansatz functions are in the multiscale space, whereas the test functions are standard finiteelement functions. The advantage over the Galerkin variant is that communication betweendifferent element correction operators is avoided. Hence, these correction operators, which arefine-scale quantities, do not need to be stored beyond the assembly of local stiffness matrixcontributions. We emphasize that in practical computations u Hm ∈ V H is first determined bysolving the linear system associated with (3.4). If required, the element correction operatorsusing u Hm can be computed to yield the full multiscale approximation u ms m . The PG-LOD (3.4) iswell-posed for sufficiently large m where no stability issues are reported in practice even for smallchoices m = 2 ,
3, cf. [9, 17]. From [9, Thm. 2] we obtain the following a priori error estimates (cid:107) A / ∇ ( u − u ms m ) (cid:107) L ( D ) + (cid:107) u − u Hm (cid:107) L ( D ) (cid:46) ( H + m d/ γ m ) (cid:107) f (cid:107) L ( D ) (3.5)for some 0 < γ < H and m . Choosing m (cid:38) | log H | , these estimates essentiallyshow that (i) u ms m converges linearly to u in the energy norm and (ii) u Hm converges linearly to u in the L ( D )-norm. In other words, u Hm is a good L -approximation to u , while the correctionoperators and thus u ms m are necessary to obtain a good H -approximation of u . In this section, we suggest an offline-online strategy for the fast computation of the left-hand sidein (3.4) for many different realizations A . For v, w ∈ V H , we denote b ( v, w ) := a ( v − C m ( A ) v, w )and observe that b ( v, w ) = (cid:88) T ∈T H b T ( v, w ) (3.6)with b T ( v, w ) := ˆ U m ( T ) A ( x )( χ T ∇ v − ∇ ( C m,T ( A ) v ))( x ) · ∇ w ( x ) dx, (3.7)where χ denotes the characteristic function. We will from now on assume that the mesh size H is an integer multiple of the periodicity length ε . This implies that b T ( · , · ) for the coefficient A ε is identical for every mesh element T ∈ T H . Hence, only b T ( · , · ) for a single T ∈ T H is requiredin order to assemble b ( · , · ) associated with A ε . Offline phase.
Fix an element T ∈ T H . Let J := { k ∈ Z d | ε ( k + Q ) ⊂ U m ( T ) } be the index setof possible defects in the patch U m ( T ) and denote by N := card J its cardinality.We introduce a bijective mapping σ : { , . . . , N } → J . Further, we set A i := (cid:40) A ε | U m ( T ) , i = 0 ,A ε | U m ( T ) + χ ε ( σ ( i )+ Q ) B ε , i = 1 , . . . , N (3.8)7 a) Random checkerboard (b) Random defect inclusions Figure 3.2.: Offline coefficients A , A , A , A (from top left to bottom right) for random checker-board of Example 2.2 (left) and random defect inclusions of Example 2.3 (right).Generated with α = 0 . β = 1 (black), ε = 2 − , H = 2 − , and m = 2.as our stored offline “basis” of coefficients. Intuitively, this means that A i is constructed from A ε by introducing a single defect. For the two examples 2.2 and 2.3 of random coefficients ofFigure 1.1 in the introduction, some corresponding A i are depicted in Figure 3.2 left and right,respectively.In the offline phase, we compute the local LOD stiffness matrix contributions b iT ( λ j , λ k ) = ˆ U m ( T ) A i ( x ) (cid:0) χ T ∇ λ j − ∇ ( C m,T ( A i ) λ j ) (cid:1) ( x ) · ∇ λ k ( x ) dx, (3.9)where { λ j } is the set of finite element basis functions spanning V H and C m,T ( A i ) denotes theelement correction operator associated with the coefficient A i . Note that the stiffness matrixcontribution for the fixed element T itself is a coarse-scale object and inexpensive to store.Additionally, we also assemble the load vector, i.e., the right-hand side of (3.4). For this assembly,we have to consider all mesh elements, but the load vector is the same for all coefficients. Online phase.
Given a sample coefficient A of the form (1.2) and (2.5), there are µ i ∈ R , i = 0 , . . . , N such that (cid:80) Ni =0 µ i = 1 and A | U m ( T ) = N (cid:88) i =0 µ i A i (3.10)for any T ∈ T H . More specifically, µ i for i = 1 , . . . , N is determined from the value of ˆ b jp fora certain j . In particular, we have µ i ∈ { , } for i = 1 , . . . , N and µ = 1 − N def where N def denotes the number of defects in the patch U m ( T ). Note that µ i depends (implicitly) on T andthat its calculation is cheap.In the online phase, we calculate the global LOD stiffness matrix as a combination of theoffline quantities as follows. With the µ i at hand, we compute the local combined LOD stiffnessmatrix contributions as ˜ b T ( λ j , λ k ) = N (cid:88) i =0 µ i b iT ( λ j , λ k ) . (3.11)8he global combined bilinear form ˜ b is defined as usual via ˜ b = (cid:80) T ∈T H ˜ b T . Effectively, the sumin (3.11) only contains N def + 1 nonzero terms. Roughly, only a fraction of p terms thus needs tobe considered each time. After the assembly of the global stiffness matrix, we compute ˜ u Hm ∈ V H as the solution of ˜ b (˜ u Hm , v H ) = F ( v H ) for all v H ∈ V H . (3.12)Note that the underlying linear system is of small dimension.More details on the (efficient) implementation of the offline-online strategy are given in Sec-tion 5.1. As already indicated by the notation, ˜ b and ˜ u Hm are only approximations to the truePG-LOD form b and solution u Hm associated with the sample coefficient A , respectively. Weanalyze the error committed by the new strategy in the next section. Remark 3.1.
Because of the periodic boundary conditions and the box-type domain, the localLOD stiffness matrix for A ε and the choice of offline coefficients are independent of T . Incase of Neumann or Dirichlet boundary conditions or more general domains, there are severalrepresentative configurations for the possible patches. One can adapt the offline phase to thesesituations by calculating and storing the matrix contributions for all possible patch configurationsand the associated offline coefficients. If the number of possible patch configurations is small,e.g., for a highly structured mesh and domain, the additional effort required may still be feasible. Remark 3.2.
The assumption of periodic A ε , B ε with H and the length of the domain beinginteger multiples of ε is important to guarantee that the offline and the samples coefficients allhave the same structure for each mesh element. Otherwise, we would need to perform the abovedescribed offline phase for all mesh elements. This, in turn, increases the computational time andthe storage costs. While the additional costs in run-time might be compensated by an effectiveonline phase if sufficiently many samples are considered, the storage may become a bottleneck.In future work, one might therefore try to reduce the number of offline coefficients per meshelement or consider adaptive online strategies, cf. Remark 4.3. Remark 3.3.
In this presentation we focus on the case when A = (cid:80) Ni =0 µ i A i . If A is notexactly a linear combination of the A i ’s one would instead need to find optimal weights { µ i } Ni =0 to minimize the error. In Section 4, we present an error estimator presented that bounds the erroralso for this case. Exactly how to do this optimization is outside the scope of this presentationand we leave it for future investigation.
4. A priori eror analysis
In this section, we discuss the well-posedness of (3.12) as well as error estimates for u − ˜ u Hm . Toaccomplish this, we start by studying the consistency error b − ˜ b . We first consider the one-dimensional case where the correction operators can be explicitly computed and then discuss thegeneralization to several dimensions. As in the previous sections, A denotes the true coefficientwith associated PG-LOD bilinear form b and the bilinear form of the offline-strategy ˜ b is definedvia (3.11).In the one-dimensional setting with I H chosen as the nodal interpolation operator, the correctorproblems automatically localize to single coarse elements, i.e., m = 0 is sufficient. Moreover, thecorrection operators can be explicitly calculated. Hence, [19] provides the following result on b T from (3.7) in d = 1: It holds for any v, w ∈ V H that b T ( v, w ) = ( A harm | T ∇ v, ∇ w ) T A harm defined as A harm | T := (cid:16) | T | ˆ T A − dx (cid:17) − . This means that b can be written as a finite element bilinear form with a modified coefficient,namely the element-wise harmonic mean. Similarly, ˜ b can be written as finite element bilinearform associated with A µ harm | T := N (cid:88) i =0 µ i A i harm | T , (4.1)where A i harm denotes the harmonic mean of A i . Let N def denote the number of defects in T for the given realization. In the following, we will write N def = θ def , T N with N the number ofpossible defect locations (in T ), cf. Section 3.3. Note that in this one-dimensional setting, we have N = H/ε . We abbreviate θ def = max T ∈T H θ def ,T . The representation of ˜ b in the one-dimensionalsetting allows for the following a priori bound. Theorem 4.1. If D is one-dimensional and I H is the nodal interpolation operator, the consis-tency error between b from (3.6) and ˜ b defined via (3.11) fulfills for any v, w ∈ V H | ( b − ˜ b )( v, w ) | ≤ βα (cid:16) β − αα (cid:17) | Q | (cid:16) εH θ def + 2 θ (cid:17) (cid:107) v (cid:107) A (cid:107) w (cid:107) A . (4.2)This theorem clearly underlines why the approach works for small defect probabilities andhence, for small θ def . We emphasize that the O ( θ def ) error contribution is multiplied by thesmall factor ε/H <
1. The proof of Theorem 4.1 is presented in Appendix A.Extending estimate (4.2) to higher dimensions is challenging because we do not have an explicitand local representation of b . In the following, we present an upper bound on the consistencyerror b − ˜ b that is computable in an a posteriori manner. We emphasize that the result does notrequire A | U m ( T ) = (cid:80) Ni =0 µ i A i . We abbreviate ¯ A = (cid:80) Ni =0 µ i A i . Theorem 4.2.
Define for any T ∈ T H E T := max v ∈ V H : v | T (cid:107) ( A / − A − / ¯ A ) χ T ∇ v − (cid:80) Ni =0 µ i ( A / − A − / A i ) ∇ ( C m,T ( A i ) v ) (cid:107) L ( U m ( T )) (cid:107) v (cid:107) A,T . (4.3) Then, for any v, w ∈ V H it holds that | (˜ b − b )( v, w ) | (cid:46) m d/ (cid:16) max T ∈T H E T (cid:17) (cid:107) v (cid:107) A (cid:107) w (cid:107) A . (4.4)In practice, E T is computed as the largest eigenvalue of an eigenvector problem of dimension2 d , which is the dimension of the coarse scale finite element space on one element T . We only needto store A i and C m,T ( A i ) on a single patch U m ( T ) and have access to the sampled coefficient A .We refer to Section 5.3 for details on these implementation aspects. Note that for N = 0, i.e., asingle “reference” coefficient A i , E T coincides with e u,T defined in [18, Lemma 3.3]. Theorem 4.2can thus be interpreted as a generalization to the case of several reference coefficients. Remark 4.3.
In the present contribution, we see E T as a computational tool to easily obtainan upper bound on the actual error without the need to compute u Hm and, in particular, the cor-rection operators C m,T ( A ) associated with A . Note that the computation of E T is not necessaryin the offline-online strategy if no error control is required. Since E T can be evaluated without C m,T ( A ), we will investigate its use to build up the offline coefficients A i or to enrich them duringthe online phase in future work. 10 roof of Theorem 4.2. Let v, w ∈ V H be arbitrary but fixed. We will show that for any T ∈ T H it holds that | ( b T − ˜ b T )( v, w ) | (cid:46) E T (cid:107) v (cid:107) A,T (cid:107) w (cid:107) A,U m ( T ) . (4.5)Let us first illustrate how this implies the assertion of the theorem: | ( b − ˜ b )( v, w ) | ≤ (cid:88) T ∈T H | ( b T − ˜ b T )( v, w ) | (cid:46) (cid:88) T ∈T H E T (cid:107) v (cid:107) A,T (cid:107) w (cid:107) A,U m ( T ) (cid:46) m d/ (cid:16) max T ∈T H E T (cid:17) (cid:107) v (cid:107) A (cid:107) w (cid:107) A . Let us now prove (4.5). We abbreviate C im,T = C m,T ( A i ) and C m,T = C m,T ( A ). We have b T ( v, w ) − ˜ b T ( v, w )= (cid:0) A ( χ T ∇ − ∇C m,T ) v, ∇ w (cid:1) U m ( T ) − N (cid:88) i =0 µ i (cid:0) A i ( χ T ∇ − ∇C im,T ) v, ∇ w (cid:1) U m ( T ) = (cid:16) ( A − ¯ A ) χ T ∇ v − N (cid:88) i =0 µ i ( A − A i ) ∇C im,T v, ∇ w (cid:17) U m ( T ) − (cid:16) A ∇ (cid:16) C m,T − N (cid:88) i =0 µ i C im,T (cid:17) v, ∇ w (cid:17) U m ( T ) ≤ E T (cid:107) v (cid:107) A,T (cid:107) w (cid:107) A,U m ( T ) + (cid:13)(cid:13)(cid:13)(cid:16) C m,T − N (cid:88) i =0 µ i C im,T (cid:17) v (cid:13)(cid:13)(cid:13) A,U m ( T ) (cid:107) w (cid:107) A,U m ( T ) . It remains to estimate the second term. We abbreviate z = C m,T v − (cid:80) Ni =0 µ i C im,T v ∈ V f ( U m ( T )).We deduce by the definition of C im,T v that (cid:107) z (cid:107) A,U m ( T ) = (cid:16) A ∇ (cid:16) C m,T v − N (cid:88) i =0 µ i C im,T v (cid:17) , ∇ z (cid:17) U m ( T ) = ( A ∇ v, ∇ z ) T − (cid:16) A N (cid:88) i =0 µ i ∇ ( C im,T v ) , ∇ z (cid:17) U m ( T ) = (( A − ¯ A ) ∇ v, ∇ z ) T − (cid:16) N (cid:88) i =0 ( A − A i ) µ i ∇C iT v, ∇ z (cid:17) U m ( T ) ≤ (cid:13)(cid:13)(cid:13) ( A / − A − / ¯ A ) χ T ∇ v − N (cid:88) i =0 µ i ( A / − A − / A i ) ∇C iT v (cid:13)(cid:13)(cid:13) L ( U m ( T )) (cid:107) z (cid:107) A,U m ( T ) ≤ E T (cid:107) v (cid:107) A,T (cid:107) z (cid:107) A,U m ( T ) , which finishes the proof.Theorems 4.1 and 4.2 provide bounds on the consistency error η := sup v ∈ V H \{ } sup w ∈ V H \{ } | ( b − ˜ b )( v, w ) |(cid:107) v (cid:107) A (cid:107) w (cid:107) A .
11f the consistency error is sufficiently small, well-posedness of (3.12) is guaranteed and we alsoobtain an error estimate as detailed in the next corollary.
Corollary 4.4.
There exist m > and η > such that, if m > m and η < η , (3.12) iswell-posed and, further, the error between the solution u of (2.2) and the solution ˜ u Hm ∈ V H of (3.12) satisfies (cid:107) u − ˜ u Hm (cid:107) L ( D ) (cid:46) ( H + m d/ γ m + η ) (cid:107) f (cid:107) L ( D ) . Proof.
We proceed similar to the proof of Theorem 4.1 in [17]. To show the well-posednessof (3.12), we prove the coercivity of ˜ b if η < η and m > m . We again abbreviate C m = C m ( A )and C m,T = C m,T ( A ). By C ∞ ,T we denote the element correction operator C m,T with m = ∞ ,i.e., where the integral on the left-hand side of (3.3) is taken over the whole domain D . Weset C ∞ = (cid:80) T ∈T H C ∞ ,T and note that a ( v − C ∞ v, C ∞ w ) = 0 for all v, w ∈ V H . Let v ∈ V H bearbitrary. We deduce˜ b ( v, v ) = b ( v, v ) + (˜ b − b )( v, v ) = a ( v − C m v, v ) + (˜ b − b )( v, v )= a ( v − C ∞ v, v ) + a ( C ∞ v − C m v, v ) + (˜ b − b )( v, v ) ≥ ( c + c m d/ γ m − η ) (cid:107) v (cid:107) A , where we used a ( v − C ∞ v, v ) = a ( v − C ∞ v, v − C ∞ v ) = (cid:107) v − C ∞ v (cid:107) A ≥ c (cid:107) I H ( v − C ∞ v ) (cid:107) A = (cid:107) v (cid:107) A and (cid:107) ( C ∞ − C m ) v (cid:107) A ≤ c m d/ γ m (cid:107) A / ∇ v (cid:107) L ( D ) in the last step. Clearly, the are m and η such that for m > m and η < η , c − c m d/ θ m − η can be bounded from below by a positive constant. This shows the coercivity of ˜ b .For the error estimate, we use the triangle inequality to get (cid:107) u − ˜ u Hm (cid:107) L ( D ) ≤ (cid:107) u − u Hm (cid:107) L ( D ) + (cid:107) u Hm − ˜ u Hm (cid:107) L ( D ) , where the first term is estimated in (3.5). By the coercivity of ˜ b we obtain for u Hm − ˜ u Hm that (cid:107) u Hm − ˜ u Hm (cid:107) A (cid:46) ˜ b ( u Hm − ˜ u Hm , u Hm − ˜ u Hm ) = ˜ b ( u Hm , u Hm − ˜ u Hm ) − F ( u Hm − ˜ u Hm )= (˜ b − b )( u Hm , u Hm − ˜ u Hm ) (cid:46) η (cid:107) f (cid:107) L ( D ) (cid:107) u Hm − ˜ u Hm (cid:107) A , where we used the stability of the PG-LOD solution u Hm in the last step. Application of Friedrich’sinequality yields the L -bound. Remark 4.5.
In the one-dimensional case, ˜ b is coercive if and only if A µ harm in (4.1) is positiveon each element. A sufficient condition – alternative to a small consistency error – is B ε ≥ A ε . In more detail, if B ε ≥
0, we can show A i harm | T ≥ A | T for i = 1 , . . . , N . Because of µ i ≥ i = 1 , . . . , N and (cid:80) Ni =0 µ i = 1, wededuce A µ harm | T ≥ α >
0. The sufficient condition B ε ≥
5. Implementation aspects
This section deals with implementation details specific to the presented offline-online strategywith a focus on the algorithm (Section 5.1), its comparison concerning run-time to the methodof [17] (Section 5.2), and the implementation of E T from Theorem 4.2 (Section 5.3). For thegeneral implementation of the (PG-)LOD we refer to [11] and [26, Ch. 7].12 .1. Algorithm for the offline-online strategy and memory consumption Algorithm 1 shows how to carry out the procedure from Section 3.3 computationally.
Algorithm 1.
Offline-online strategy input: Problem data A ε , B ε , Q , f (cf. Section 2) Pick m Fix T ∈ T H (cid:46) start offline phase Precompute and save offline coefficients { A , A , . . . A N } according to (3.8) for i = 0 , . . . N do Precompute C m,T ( A i ) λ j for all j (discard at end of iteration) Precompute and save b iT ( λ j , λ k ) for all j and k end for Precompute and save F ( λ j ) for all j (cid:46) end offline phase for all sample coefficients A do (cid:46) start online phase for all T ∈ T H do Compute µ i such that A | U m ( T ) = (cid:80) Ni =0 µ i A i Compute and save ˜ b T ( λ j , λ k ) = (cid:80) Ni =0 µ i b iT ( λ j , λ k ) for all j and k end for Assemble stiffness matrix ˜ K kj := (cid:80) T ∈T H ˜ b T ( λ j , λ k ) Solve for ˜ u Hm according to (3.12) end for It starts with setting up the offline coefficients. Due to the periodicity and the weakly randomstructure, only the fine-scale representations of A ε and B ε on a single coarse element T need tobe stored with a cost of order ( H/h ) d . b iT is computed for all offline coefficients, but only for asingle coarse element. The storage cost is of order N m d where the number of offline coefficients N is of the order ( mH/ε ) d . We especially emphasize that, for moderate m and not too coarse H , we have N < O ( ε − d ) = card I with the index set I from (2.5). This means that there areless possible defects in the element patch U m ( T ) than in D . The pre-computations for the offlinecoefficients can be executed in parallel. Finally, we also assemble and store the load vector inthe offline phase with a cost of order H − d .In the online phase, we perform a loop over the sample coefficients (i.e., Monte-Carlo-typesampling) which again can be executed in parallel. For each coefficient, there is a loop overthe elements which is also parallelizable. For each element, we extract the µ i forming therepresentation of A in terms of the offline coefficients. Due to the representation (2.5), thisdoes not require extensive computations. We emphasize that many µ i are identically zero forrare perturbations so that the sum for ˜ b T is evaluated cheaply. Finally, the coarse-scale linearsystem with ˜ b is assembled and solved. The assembly of the stiffness matrix involves a reductionover T , but only coarse-scale quantities like ˜ b T of amount m d H − d are needed between differentelements, cf. [17]. From the coarse-scale solution ˜ u Hm for each sample coefficient, we can of coursecompute quantities of interest or agglomerate statistical information. Based on Algorithm 1, let us briefly comment on the run-time complexity in comparison to thestandard LOD and the LOD with local updates according to [17]. We consider the time taken forthe stiffness matrix assembly for M samp sample coefficients and consider completely sequentialversions of all methods. In the following, t stiff measures the time for the assembly of a localLOD stiffness matrix contribution for a given coefficient according to (3.7). We assume that this13ime does not depend on the given coefficient or the coarse element T considered. Further, n H denotes the number of elements in T H which is of the order H − d .For the standard LOD, the global LOD matrix is newly computed for each sample. Hence, thetotal time for LOD matrix assemblies amounts to t s tot := M samp n H t stiff . In the LOD with localupdates of [17], the LOD stiffness matrix for A = A ε is calculated on a single element. For eachsample coefficient A , an error indicator is then evaluated for each element T which requires atime of n H t ind . For the fraction p recomp of elements where the error indicator is the largest, theLOD stiffness matrix is computed anew based on A . Overall, the LOD with local updates takesa total matrix assembly time of t u tot := t stiff + M samp ( n H t indic + p recomp n H t stiff ) . In the offline-online strategy, we first compute LOD stiffness matrices for ( N + 1) coefficients,yielding a run-time of ( N + 1) t stiff in the offline phase. We recall that N is of the order ( mH/ε ) d .In the online phase, we denote by t comb the time to combine the LOD stiffness matrices. Hence,we obtain the total matrix assembly time for the offline-online strategy as t o tot := ( N + 1) t stiff + M samp n H t comb . We observe that the offline-online strategy easily outperforms the standard LOD, i.e., t o tot < t s tot because we can expect t comb (cid:28) t stiff : Forming the linear combination in (3.12) is much fasterthan computing the correction operator C m,T ( A ). The main goal thus is to outperform the LODwith local updates, i.e., to achieve t o tot < t u tot . This requires t comb < t indic + p recomp t stiff , whichis easily achievable in practice, and we can even hope for t comb < t indic . We then deduce that t o tot < t u tot if and only if M samp > Nt stiff n H ( t indic + p recomp t stiff − t comb ) . With the previous comments onthe relation of t comb , t stiff and t ind , the offline-online strategy will be more efficient than theLOD with local updates already for moderate sample sizes. This holds especially in the regime H ≈ √ ε where N/n H is small. E T We finally give some details on the implementation of E T from Theorem 4.2, where we onlyconsider the case A = ¯ A for simplicity. When computing E T we set up an eigenvalue problem.We denote by λ j the local basis functions of V H on T . There are d + 1 basis functions forsimplicial meshes and 2 d on quadrilateral meshes. Now, E T is the square root of the largesteigenvalue to the generalized eigenvalue problem S v = ν B v, where S jk = (cid:16) N (cid:88) i =0 µ i ( A / − A − / A i ) ∇ ( C m,T ( A i ) λ k ) , N (cid:88) i =0 µ i ( A / − A − / A i ) ∇ ( C m,T ( A i ) λ j ) (cid:17) U m ( T ) B jk = ( A ∇ λ k , ∇ λ j ) T . To assemble S , we need to store ∇C m,T ( A i ) λ j for i = 0 , . . . , N and j = 1 , . . . , d (since (cid:80) λ j =1) on each (fine) element of the patch surrounding T (in case of a simplicial mesh). Since thesegradients are piecewise constant on the fine mesh, this amounts to ( N + 1) d vectors of length d per fine element. Note that the number of fine elements in a patch is of the order ( mH/h ) d .Hence, additional to the cheap storage of A i (see Section 5.1), quantities of order N d ( mH/h ) d have to be stored for computing E T . We emphasize that ∇C m,T ( A i ) λ j only needs to be stored fora single, fixed element T due to periodicity. For each A = (cid:80) Ni =0 µ i A i we then need to computethe component-wise products between ( A / − A − / A i ) µ i and the stored ∇C m,T ( A i ) λ j , sumover i and finally compute the integrals for the combinations of 1 ≤ j, k ≤ d + 1. Here, we canexploit symmetry and the fact that many µ i are zero. We emphasize that the computation of B is cheaper. 14 emark 5.1. We note once more that E T is not required in the offline-online strategy, butserves as error control. Hence, instead of evaluating E T for each sample coefficient, one couldalso try to estimate the distribution of E T by sampling (via sampling µ i ). This could be donefor a single element T in the (extended) offline phase so that C m,T ( A i ) can be discarded for theonline phase as before.
6. Numerical experiments
In this section, we present extensive numerical examples in one and two dimensions on theunit cell D = [0 , d . We choose f = 8 π sin(2 πx ) in the one-dimensional experiment and f = 8 π sin(2 πx ) cos(2 πx ) in two dimensions. Our code is based on gridlod [16] and ispublicly available at https://github.com/BarbaraV/gridlod-random-perturbations .We mainly focus on the relative L -error (cid:107) u Hm − ˜ u Hm (cid:107) L ( D ) (cid:107) u Hm (cid:107) L ( D ) (6.1)between the PG-LOD solution u Hm for the given sample coefficient and the solution ˜ u Hm of ouroffline-online strategy. We additionally take the root mean square error over M samp samples. We set h = ε = 2 − and randomly assign to each interval of length ε the value α = 0 . − p or the value β = 1 . p . We compute the maximal error betweenthe element-wise harmonic means A harm | T and A µ harm | T , cf. Section 4, as well as the relative L ( D )-errors for the solutions (with m = 0). The root mean square errors over 500 samples independence on the probability p and on the mesh size H are depicted in Figure 6.1. On the left,we see that all curves for the harmonic means show an overall quadratic behavior as expectedfrom Theorem 4.1. Moreover, we also see the predicted increase of the constant with decreasing H . Figure 6.1 (right) illustrate the quadratic dependence of the root mean square errors in thesolution on p . Here, however, the curves for different H lie closer together. Moreover, we observethat the relative L ( D )-errors (slightly) decrease with decreasing H . Overall, we clearly see thatup to probabilities of p = 0 .
2, our new method produces root mean square errors (in the solution)of less than 3% which is acceptable in many applications.
We now consider a two-dimensional random checkerboard coefficient, i.e., we set d = 2 andconsider Example 2.2 with ε = 2 − , α = 0 . β = 1 .
0. For the LOD we set h = 2 − , sothat all fine-scale details are resolved, as well as H = 2 − and m = 4. The choice of H and m ensures that u Hm of the standard LOD is a good reference solution. All the following resultsare obtained with M samp = 350. First, we investigate the behavior of the root mean squarerelative L ( D )-error of our new approach in dependence on the probability p . We observe inFigure 6.2 (left) again a quadratic growth as in the one-dimensional setting. Note that the rootmean square error stays below roughly 3% for probabilities up to p = 0 . A ε – which is deterministic – gives a rather poor approximation as expected, see Figure 6.2,right. More precisely, the root mean square error grows linearly in p and is already about 10%for p = 0 .
01. This also implies that the LOD with local updates from [17] needs an increasingfraction of basis updates. For instance, for p = 0 .
1, we would require about 60% updates in15igure 6.1.: Root mean square errors for harmonic means in L ∞ ( D )-norm (left) and root meansquare relative L ( D )-errors of the solution (right) versus probability for the one-dimensional random checkerboard with ε = 2 − and varying H .Figure 6.2.: Root mean square relative L ( D )-errors versus probability for the two-dimensionalrandom checkerboard. Left: offline-online strategy, right: LOD solution for deter-ministic coefficient A ε .this example to attain an accuracy comparable to the offline-online strategy. Put differently,for p = 0 .
1, the LOD with 15% updates (which is a reasonable amount and considered for thetimings below) leads to a root mean square error of about 17%, in contrast to the reported 3%for the offline-online strategy.In the spirit of Section 5.2, let us briefly comment on the run-time of our offline-online-strategy.We consider the same setting as before and fix p = 0 .
1. Timings were taken without theparallelization possibilities discussed in Section 5.1 on a standard desktop computer (Intel i5-7500 core, frequency 3.40 GHz, Ubuntu 20.04). With the proposed strategy, the offline phasetakes about 186 seconds and the assembly of the global stiffness matrix takes less than 3 secondsfor a single sample coefficient in the online phase (averaged over 350 samples). In contrast,computing a new LOD stiffness matrix completely for each coefficient takes about 145 secondsper sample (averaged over 350 samples). This clearly shows that a standard LOD becomesextremely costly in many Monte Carlo settings. Computing the LOD stiffness matrix on a singleelement for A ε takes about 0 . ε = 2 − and p = 0 .
1. Top row: Changes in value with ˜ β ∈ { ., . , } (from left to right).Bottom row from left to right: fill , shift , and Lshape .as matrix assembly time – including evaluation of the corrector and local re-computations – for asingle sample (averaged over 350 samples). This indicates that the offline-online strategy is moreattractive than local updates already for moderate sample sizes – in the considered example,after about 10 samples.
We now consider the two-dimensional version of Example 2.3 with ε = 2 − , α = 1 and β = 10.For the PG-LOD, we select h = 2 − , H = 2 − and m = 3 again guaranteeing that all inclusionsare resolved by the fine mesh and that the PG-LOD solution u Hm serves as a good reference.We consider the root mean square relative L ( D )-errors for our offline-online strategy over 350samples and for defect probabilities p ∈ { . , . , . , . } . Here we focus on the influence ofthe following defect possibilities on the errors: • A defect inclusion has the new value ˜ β ∈ { , . , } where we emphasize that ˜ β = 1 meansthat the inclusion vanishes; • A defect inclusion takes the value β in the whole ε -cell (called fill ); • A defect inclusion is positioned at (scaled and shifted versions of) [0 . , (called shift );17igure 6.4.: Root mean square relative L ( D )-error versus probability for the two-dimensionalperiodic inclusions with defects. Comparison of different defect values (left) andgeometry changes (right). • A defect inclusion has a different shape of (scaled and shifted versions of) [0 . , . \ [0 . , . (called Lshape ).The different considered possibilities for p = 0 . ε = 2 − to better see the (fine) inclusions. Note that the model shift does not only shift theinclusion, but even changes its size.The root mean square relative L ( D )-errors of our offline-online strategy depending on p aredepicted in Figure 6.4. On the left, we see the influence of the value taken in a defect. While allroot mean square errors are very small (below 0 . β , i.e., the value in the defect,the larger the error. The dependency of the error on the contrast between α , β , and ˜ β is alsoindicated by our theoretical findings.Figure 6.4 (right) shows the influence of the different geometry changes in the defects on theroot mean square errors. We emphasize that all errors of the new strategy stay below 5% upto p = 0 .
15. Concerning the differences between the geometry changes, we clearly observe thatthe model fill is the most difficult and produces the largest errors. This holds true not only incomparison to shift and
Lshape depicted in Figure 6.4 right, but also in comparison to value depicted in Figure 6.4 (left). E T The aim of our final numerical experiment is to illustrate the relation between the error u Hm − ˜ u Hm locally on an element patch U m ( T ) and the indicator E T . Note that this is not exactly coveredby the theoretical findings of Theorem 4.2. We approach the question by studying the randomcheckerboard coefficient of Example 2.2 in d = 2 with α = 0 . β = 1 .
0. We set H = 1 / m = 2, ε = 1 /
20 and h = 1 /
40 so that D = [0 , represents the element patch U m ( T ). For eachsample coefficient, we compute u Hm and ˜ u Hm as in the previous experiments. In this experiment,we consider both the root mean square absolute and relative L ( D )-error. Additionally, wecompute E T for the element T in the middle of D . Figure 6.5 depicts the root mean squares –sampled over 500 realizations – of the errors and the indicator for different values of p . We seethat the absolute L ( D )-error and E T show a qualitatively and quantitatively similar behaviorwhich underlines the validity of E T as an indicator for the (local) error. The relative L ( D )-erroris of a different magnitude since we divide by the norm of u Hm , cf. (6.1). Amplifying the relative18igure 6.5.: Root mean square of relative and absolute L ( D )-errors and indicator E T versusprobability for the two-dimensional random checkerboard.error by a factor of 4, we observe that there seems to be a good qualitative agreement betweenthe relative error and the indicator (Figure 6.5, dashed line). In other words, the ratio between E T and the absolute as well as the relative L ( D )-error seems to be almost constant for varyingprobabilities p . Conclusion
We presented an offline-online strategy based on the Localized Orthogonal Decomposition (LOD)to compute coarse-scale solutions of elliptic equations with local random perturbations in a MonteCarlo setting. Exploiting the periodic structure of the underlying deterministic coefficient, LODstiffness matrices on a single reference element for representative perturbations are computedin an offline phase. Those are efficiently combined to yield the LOD stiffness matrix for eachsample coefficient in the online phase. We showed error estimates in the one- as well as the higherdimensional case where we derived an indicator. Theoretical algorithmic considerations as wellas numerical experiments underlined the promising performance of the offline-online strategy.Overall, the present contribution provides interesting first results on an efficient computationalmultiscale approach for PDEs with random coefficients. Future research concerns the combi-nation with Monte Carlo-type approaches, the extension of the construction to more generalrandom coefficients, and the improvement of the error indicator to reduce storage of fine-scalequantities.
References [1] A. Abdulle and P. Henning. A reduced basis localized orthogonal decomposition.
J. Comput.Phys. , 295:379–401, 2015.[2] A. Anantharaman and C. Le Bris. A numerical approach related to defect-type theories for19ome weakly random problems in homogenization.
Multiscale Model. Simul. , 9(2):513–544,2011.[3] A. Anantharaman and C. Le Bris. Elements of mathematical foundations for numeri-cal approaches for weakly random homogenization problems.
Commun. Comput. Phys. ,11(4):1103–1143, 2012.[4] I. Babuˇska, F. Nobile, and R. Tempone. A stochastic collocation method for elliptic partialdifferential equations with random input data.
SIAM J. Numer. Anal. , 45(3):1005–1034,2007.[5] A. Barth, C. Schwab, and N. Zollinger. Multi-level Monte Carlo finite element method forelliptic PDEs with stochastic coefficients.
Numer. Math. , 119(1):123–161, 2011.[6] X. Blanc, C. Le Bris, and F. Legoll. Some variance reduction methods for numerical stochas-tic homogenization.
Philos. Trans. Roy. Soc. A , 374(2066):20150168, 2016.[7] K. A. Cliffe, M. B. Giles, R. Scheichl, and A. L. Teckentrup. Multilevel Monte Carlo methodsand applications to elliptic PDEs with random coefficients.
Comput. Vis. Sci. , 14(1):3–15,2011.[8] Y. Efendiev, C. Kronsbein, and F. Legoll. Multilevel Monte Carlo approaches for numericalhomogenization.
Multiscale Model. Simul. , 13(4):1107–1135, 2015.[9] D. Elfverson, V. Ginting, and P. Henning. On multiscale methods in Petrov-Galerkin for-mulation.
Numer. Math. , 131(4):643–682, 2015.[10] D. Elfverson, F. Hellman, and A. M˚alqvist. A multilevel Monte Carlo method for computingfailure probabilities.
SIAM/ASA J. Uncertain. Quantif. , 4(1):312–330, 2016.[11] C. Engwer, P. Henning, A. M˚alqvist, and D. Peterseim. Efficient implementation of thelocalized orthogonal decomposition method.
Comput. Methods Appl. Mech. Engrg. , 350:123–153, 2019.[12] M. Feischl and D. Peterseim. Sparse Compression of Expected Solution Operators.
SIAMJ. Numer. Anal. , 58(6):3144–3164, 2020.[13] J. Fischer, D. Gallistl, and D. Peterseim. A priori error analysis of a numerical stochastichomogenization method. arXiv preprint 1912.11646, 2019. Accepted for publication in
SIAM J. Numer. Anal. [14] D. Gallistl and D. Peterseim. Numerical stochastic homogenization by quasilocal effectivediffusion tensors.
Commun. Math. Sci. , 17(3):637–651, 2019.[15] M. D. Gunzburger, C. G. Webster, and G. Zhang. Stochastic finite element methods forpartial differential equations with random input data.
Acta Numer. , 23:521–650, 2014.[16] F. Hellman and T. Keil. Gridlod. https://github.com/fredrikhellman/gridlod , 2019.GitHub repository.[17] F. Hellman, T. Keil, and A. M˚alqvist. Numerical upscaling of perturbed diffusion problems.
SIAM J. Sci. Comput. , 42(4):A2014–A2036, 2020.[18] F. Hellman and A. M˚alqvist. Numerical homogenization of elliptic PDEs with similar coef-ficients.
Multiscale Model. Simul. , 17(2):650–674, 2019.2019] P. Hennig, R. Maier, D. Peterseim, D. Schillinger, B. Verf¨urth, and M. K¨astner. Adiffuse modeling approach for embedded interfaces in linear elasticity.
GAMM-Mitt. ,43(1):e202000001, 2020.[20] P. Henning and D. Peterseim. Oversampling for the multiscale finite element method.
Multiscale Model. Simul. , 11(4):1149–1175, 2013.[21] C. Le Bris. Homogenization theory and multiscale numerical approaches for disorderedmedia: some recent contributions. In
Congr`es SMAI 2013 , volume 45 of
ESAIM Proc.Surveys , pages 18–31. EDP Sci., Les Ulis, 2014.[22] C. Le Bris, F. Legoll, and F. Thomines. Multiscale finite element approach for “weakly”random problems and related issues.
ESAIM Math. Model. Numer. Anal. , 48(3):815–858,2014.[23] C. Le Bris and F. Thomines. A reduced basis approach for some weakly stochastic multiscaleproblems.
Chin. Ann. Math. Ser. B , 33(5):657–672, 2012.[24] G. J. Lord, C. E. Powell, and T. Shardlow.
An introduction to computational stochasticPDEs . Cambridge Texts in Applied Mathematics. Cambridge University Press, New York,2014.[25] A. M˚alqvist and D. Peterseim. Localization of elliptic multiscale problems.
Math. Comp. ,83(290):2583–2603, 2014.[26] A. M˚alqvist and D. Peterseim.
Numerical Homogenization by Localized Orthogonal De-composition . SIAM Spotlights. Society for Industrial and Applied Mathematics (SIAM),Philadelphia, PA, 2020.[27] M. Ohlberger and B. Verf¨urth. Localized Orthogonal Decomposition for two-scaleHelmholtz-type problems.
AIMS Mathematics , 2(3):458–478, 2017.[28] N. Ou, G. Lin, and L. Jiang. A Low-Rank Approximated Multiscale Method for Pdes WithRandom Coefficients.
Multiscale Model. Simul. , 18(4):1595–1620, 2020.[29] A. L. Teckentrup, R. Scheichl, M. B. Giles, and E. Ullmann. Further analysis of multilevelMonte Carlo methods for elliptic PDEs with random coefficients.
Numer. Math. , 125(3):569–600, 2013.[30] Z. Zhang, M. Ci, and T. Y. Hou. A multiscale data-driven stochastic method for ellipticPDEs with random coefficients.
Multiscale Model. Simul. , 13(1):173–204, 2015.
A. Proof of Theorem 4.1
For the proof of Theorem 4.1, we use the notation from Section 4 and introduce some furtherabbreviations. We recall the coefficient A µ harm | T associated with ˜ b T for given µ i with (cid:80) Ni =0 µ i = 1,see (4.1). Further, we introduce A := ˆ ε A ε dx (A.1)and point out that we have ˆ T A ε dx = N A
21y the periodicity of A ε . Similarly, we denote A def = ˆ ε ( k + Q ) A ε + B ε − A ε dx (A.2)and emphasize that A def is independent of the choice of k ∈ I with the index set I from (2.5). Proof of Theorem 4.1.
We estimate the consistency error as | ( b − ˜ b )( v, w ) | ≤ (cid:88) T ∈T H | ( b T − ˜ b T )( v, w ) |≤ (cid:88) T ∈T H (cid:12)(cid:12) A harm | T − A µ harm | T (cid:12)(cid:12) (cid:107)∇ v (cid:107) L ( T ) (cid:107)∇ w (cid:107) L ( T ) ≤ α − (cid:0) max T ∈T H (cid:12)(cid:12) A harm | T − A µ harm | T (cid:12)(cid:12)(cid:1) (cid:107) v (cid:107) A (cid:107) w (cid:107) A . Hence, it suffices to estimate | A harm | T − A µ harm | T | for any T ∈ T H .With the preliminary observations and notation from above, we have A harm | T = | T | N A + N def A def ,A | T = | T | N A , and A i harm | T = | T | N A + A def for i = 1 , . . . , N. Recall that for i = 1 , . . . , N we have µ i ∈ { , } and µ = 1 − N def . Hence, we deduce A µ harm | T = (cid:88) µ i A i harm | T = (1 − N def ) | T | N A + N def | T | N A + A def = N A | T | + (1 − N def ) A def | T | N A ( N A + A def )= | T | N A − N def A def | T | N A ( N A + A def ) . (A.3)Combining (A.3) with the expression for A harm | T , inserting N def = θ def ,T N , and performing aTaylor expansion around θ = 0, we obtain | A harm | T − A µ harm | T | = | T | (cid:12)(cid:12)(cid:12) N A + θ def , T N A def − N A + θ def ,T A def A ( N A + A def ) (cid:12)(cid:12)(cid:12) ≤ | T | (cid:16) θ def ,T (cid:12)(cid:12)(cid:12) A def A ( N A + A def ) − A def N A (cid:12)(cid:12)(cid:12) + θ ,T (cid:12)(cid:12)(cid:12) A N ( A + ηA def ) (cid:12)(cid:12)(cid:12)(cid:17) for some η ∈ [0 , θ def , T ].Before we estimate the first- and second-order term in the expansion separately, we bound A and A def using (2.3)–(2.4). We obtain εβ ≤ A ≤ εα and A def ≤ ε | Q | (cid:16) α − β (cid:17) as well as (writing Q = [ q , q ]) N A + A def = N ˆ ε A ε dx + ˆ εq εq A ε + B ε − A ε dx
22 ( N − ˆ ε A ε dx + ˆ εq A ε dx + ˆ εεq A ε dx + ˆ εq εq A ε + B ε dx ≥ β (( N − ε + εq + ε (1 − q ) + ε ( q − q )) = N εβ .
For the first-order term in the Taylor expansion we deduce (cid:12)(cid:12)(cid:12) A def A ( N A + A def ) − A def N A (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12) A def ( N A − ( N A + A def )) N A ( N A + A def ) (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12) A N A ( N A + A def ) (cid:12)(cid:12)(cid:12) ≤ ε | Q | (cid:16) α − β (cid:17) β ε N = β | Q | ε N (cid:16) α − β (cid:17) . Similarly, the second-order term in the Taylor expansion can be estimated as (cid:12)(cid:12)(cid:12) A N ( A + ηA def ) (cid:12)(cid:12)(cid:12) ≤ ε | Q | N (cid:16) α − β (cid:17) β ε (1 + | Q | (1 − η )) ≤ ε | Q | β N ε (cid:16) α − β (cid:17) . Finally, we obtain with | T | = H and N = H/ε that | A harm | T − A µ harm | T | ≤ | Q | β (cid:16) α − β (cid:17) εH θ def ,T + 2 | Q | β (cid:16) α − β (cid:17) θ ,T ..