a posteriori analysis and adaptive error control for multiscale operator ...
A POSTERIORI ANALYSIS AND ADAPTIVE ERROR CONTROL FOR MULTISCALE OPERATOR DECOMPOSITION SOLUTION OF ELLIPTIC SYSTEMS I: TRIANGULAR SYSTEMS V. CAREY∗ , D. ESTEP† , AND S. TAVENER‡ Abstract. In this paper, we perform an a posteriori error analysis of a multiscale operator decomposition finite element method for the solution of a system of coupled elliptic problems. The goal is to compute accurate error estimates that account for the effects arising from multiscale discretization via operator decomposition. Our approach to error estimation is based on a well known a posteriori analysis involving variational analysis, residuals and the generalized Green’s function. Our method utilizes adjoint problems to deal with several new features arising from the multiscale operator decomposition. In Part I of this paper, we focus on the propagation of errors arising from the solution of one component to another and the transfer of information between different representations of solution components. We also devise an adaptive discretization strategy based on the error estimates that specifically controls the effects arising from operator decomposition. In Part II of this paper, we address issues related to the iterative solution of a fully coupled nonlinear system. Key words. a posteriori error analysis, adjoint problem, elliptic system, generalized Green’s function, goal oriented error estimates, multiscale methods, operator decomposition, projection error AMS subject classifications. 65N15, 65N30, 65N50
1. Introduction. Multiscale operator decomposition is a widely used technique for solving multiphysics, multiscale problems. The general approach is to decompose the multiphysics problem into components involving simpler physics over a relatively limited range of scales, and then to seek the solution of the entire system through some sort of iterative procedure involving solutions of the individual components. This approach is appealing because there is generally a good understanding of how to solve a broad spectrum of single physics problems accurately and efficiently, and because it provides an alternative to accommodating multiple scales in one discretization. However, multiscale operator decomposition presents an entirely new set of accuracy and stability issues, some of which are obvious and some subtle, and all of which are difficult to correct. We motivate multiscale operator decomposition for elliptic systems by considering a model of a thermal actuator. A thermal actuator is a MEMS (microelectronic mechanical switch) device. A contact rests on thin braces composed of a conducting material. When a current is passed through the braces, they heat up and consequently expand to close the contact. The system is modelled by a system of three coupled equations, each representing a distinct physical process. They are an electrostatic ∗ Department of Mathematics, Colorado State University, Fort Collins, CO 80523 ([email protected]). V. Carey’s work is supported in part by the Department of Energy (DE-FG02-04ER25620) † Department of Mathematics and Department of Statistics, Colorado State University, Fort Collins, CO 80523 ([email protected]). D. Estep’s work is supported in part by the Department of Energy (DE-FG02-04ER25620, DE-FG02-05ER25699, DE-FC02-07ER54909), the National Aeronautics and Space Administration (NNG04GH63G), the National Science Foundation (DMS0107832, DMS-0715135, DGE-0221595003, MSPA-CSE-0434354, ECCS-0700559), Idaho National Laboratory (00069249), and the Sandia Corporation (PO299784) ‡ Department of Mathematics, Colorado State University, Fort Collins, CO 80523 ([email protected]). S. Tavener’s work is supported in part by the Department of Energy (DE-FG02-04ER25620)
1
2
V. CAREY, D. ESTEP, AND S. TAVENER Vsig
Vswitch
Fig. 1.1: Sketch of a thermal actuator.
current equation ∇ · (σ∇V ) = 0,
(1.1)
governing potential V (where current J = −σ∇V ), a steady-state energy equation ∇ · (κ(T )∇T ) = σ(∇V · ∇V ),
(1.2)
governing temperature T , and a linear elasticity equation giving the steady-state displacement d, ¡ ¢ ¡ ¢ ∇ · λ tr(E)I + 2µE − β(T − Tref )I = 0, E = ∇d + ∇d> /2. (1.3) Using multiscale operator decomposition, the complete system (1.1-1.3) is decomposed into three components, each of which is solved with a code specialized to the particular type of physics. Notice that the electric potential V can be calculated independently of T and d. The temperature T can be calculated once the electric potential V is known, while the calculation of displacement d requires prior knowledge of T , and therefore of V . In general, we can write a coupled elliptic system on a domain Ω in the form L (x, u1 , Du1 , . . . , un , Dun ) = 0 1 .. , x ∈ Ω. (1.4) . L (x, u , Du , . . . , u , Du ) = 0 n 1 1 n n A natural form of operator decomposition is to split the global multi-physics problem into n “single-physics” components that are solved individually. In general, the solution of each component requires knowledge of the solutions of all the other components; the full problem requires some form of iteration to obtain the solution. It is possible to impose conditions on the system, the components, and the coupling that allow for an a priori convergence analysis. However, operator decomposition is problematic in practice because it is very difficult to verify such conditions
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
3
and often impractical to satisfy them. Indeed, numerical solutions obtained via operator decomposition are affected significantly by the specific choice of decomposition. In this paper, we perform an a posteriori error analysis of a multiscale operator decomposition finite element method for the solution of a system of coupled elliptic problems. The components of the problem are solved in sequence using independent discretizations. The goal is to compute accurate computational error estimates that specifically account for the effects arising from operator decomposition. We also devise an adaptive discretization strategy based on the error estimates that controls the effects arising from multiscale operator decomposition. The a posteriori analysis in this paper is based on a well known approach involving variational analysis, residuals and the generalized Green’s function solving an adjoint problem [1, 2, 5–7, 9, 12]. However, we modify this approach to accommodate several new features arising from the operator decomposition. Three important issues addressed here are: (1) errors in the solution of each component propagate into the solutions of the other components; (2) transferring information between different discretization representations potentially introduces new error; and (3) the adjoint operators associated with the fully coupled system and an operator decomposition version are not generally equal. In addition, the analysis stays within the “single physics paradigm” by only requiring the solution of adjoint problems associated with the individual components. These issues are characteristic of a broad range of operator decomposition discretizations, e.g. [10,13], and generally require extensions to the usual a posteriori analysis techniques. In this paper, we focus attention on analyzing the effects of transferring information between components, which is necessitated by operator decomposition. In order to do so, we consider a “triangular” or one-way coupled system, L1 (x, u1 , Du1 ) = 0 L (x, u1 , Du1 , u2 , Du2 ) = 0 2 L3 (x, u1 , Du1 , u2 , Du2 , u3 , Du3 ) = 0 , x ∈ Ω. (1.5) .. . Ln (x, u1 , Du1 , u2 , Du2 , u3 , Du3 , . . . , un , Dun ) = 0 This system can be solved by a finite sequence of component solutions by considering the n problems for L1 , L2 , · · · , Ln sequentially. Such systems are important in practice, e.g., the thermal actuator (1.1)–(1.3) has this form. In Part II [3], we consider additional sources of error arising from the iterative procedure required when solving a fully-coupled” system via operator decomposition. We capture the essential features of (1.5) in a two component “one-way” coupled system of the form x ∈ Ω, −∇ · a1 ∇u1 + b1 · ∇u1 + c1 u1 = f1 (x), (1.6) −∇ · a2 ∇u2 + b2 · ∇u2 + c2 u2 = f2 (x, u1 , Du1 ), x ∈ Ω, u1 = u2 = 0, x ∈ ∂Ω, where ai , bi , ci , fi are smooth functions on a bounded domain Ω in RN with boundary ∂Ω, and the coupling occurs through f2 . We later generalize to coupling through the coefficients of the elliptic operator for u2 . In Section 2, we illustrate the main idea by applying the analysis to a linear algebraic system. We perform the transfer error analysis in Section 3 and present
4
V. CAREY, D. ESTEP, AND S. TAVENER
computational examples when the corresponding discretizations are “related” in the sense that either both computational meshes are identical, or one mesh is generated by a sequence of mesh refinements on the other mesh. In Section 4, we consider the effect of using distinct discretizations for the two components and analyze the additional errors cause by using projections between the components. Additionally, we discuss the use of Monte Carlo integration to estimate these projection errors. We present the full adaptive algorithm in section 5 which we illustrate with several numerical examples. 2. A linear algebra example. We introduce the notation and ideas in the context of a lower triangular linear system of equations. Let U be an approximate solution of the linear system Au = b. We wish to compute a quantity of interest given by a linear functional (ψ, u). The error, e = u − U is not computable, but we can compute the residual R = b − AU = Ae. Using the solution φ of the corresponding adjoint equation A> φ = ψ, the error representation for a linear functional of the solution is ¡ ¢ ¡ ¢ ¡ ¢ ¡ ¢ ¡ ¢ ¡ ¢ ψ, u − ψ, U = ψ, e = A> φ, e = φ, Ae = φ, R . Now consider the triangular system µ ¶µ ¶ µ ¶ A11 0 u1 b Au = = 1 =b, A21 A22 u2 b2
(2.1)
with approximate solution µ ¶ µ ¶ U1 u1 U= ≈ = u. U2 u2 We estimate the error in a quantity of interest in u2 only, given by the linear functional µ ¶ ¢ ¡ (1) ¢ ¡ (1) 0 . ψ , u = ψ2 , u2 where ψ = (1) ψ2 We employ the superscript (1) since we later pose additional auxiliary adjoint problems. Clearly, estimates on linear functionals of u1 are independent of u2 . The lower triangular structure of A yields A11 u1 = b1 , A22 u2 = b2 − A21 u1 , and the corresponding residuals are R1 = b1 − A11 U1 , R2 = (b2 − A21 U1 ) − A22 U2 . The residual R2 depends upon the solution of the first component, and any attempt to decrease this residual requires a consideration of the accuracy of U1 . The adjoint problem to (2.1) is µ > ¶ Ã (1) ! µ ¶ 0 A11 A> φ1 21 = (1) , (1) 0 A> ψ2 φ2 22
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
and the resulting error representation is ¡ (1) ¡ (1) ¢ ¡ ¢ (1) ψ , e) = ψ2 , e2 = A> 22 φ2 , e2 ¢ ¢ ¡ (1) ¡ (1) = φ2 , A22 u2 − φ2 , A22 U2 ¡ (1) ¢ ¡ (1) ¢ = φ2 , b2 − A21 u1 − φ2 , A22 U2 ¡ (1) ¢ ¡ (1) ¢ = φ2 , b2 − A21 U1 − A22 U2 + φ2 , A21 e1 ¡ (1) ¢ ¡ (1) ¢ = φ2 , R2 + φ2 , A21 e1 .
5
(2.2)
(1)
The first term of the error representation requires only U2 and φ2 . Since the adjoint system is upper triangular and ³ ´−1 (1) (1) φ 2 = A> ψ2 22 ¡ (1) ¢ is independent of the first component, the calculation of φ2 , R2 remains within ¡ (1) ¢ the “single physics paradigm”. The second term φ2 , A21 e1 represents the effect of errors in U1 on the solution U2 . At first glance this term is uncomputable, but we note that it is a linear functional of e1 since ¢ ¢ ¡ ¡ (1) (1) φ2 , A21 e1 = A> 21 φ2 , e1 . We therefore form the adjoint problem for the transfer error, ¶ Ã (2) ! µ (2) ¶ µ > (1) ¶ µ > φ1 A11 A> ψ1 A21 φ2 21 = = . (2) 0 A> 0 0 φ 22 2 (2)
The upper triangular block structure of A> immediately yields φ2 = 0. As noted (2) (2) earlier, error estimates of u1 should be independent of u2 . Thus, A> = 11 φ1 = ψ1 > (1) (2) A21 φ2 , so that once again we can solve for φ in the “single physics paradigm”. Given φ(2) we obtain the secondary error representation ¢ ¡ (2) ¢ ¡ (2) ¢ ¡ > (1) ¢ ¡ > (2) ¢ ¡ (2) (2.3) ψ , e = ψ1 , e1 = A21 φ2 , e1 = A11 φ1 , e1 = φ1 , R1 . Combining the first term of (2.2) with (2.3) yields the complete error representation ¡ (1) ¢ ¡ (1) ¢ ¡ (2) ¢ ψ , e = φ2 , R2 + φ1 , R1 (2.4) which is a sum of the inner products of “single physics” residuals and adjoint solutions computed using the “single physics” paradigm. 3. Analysis of the discretization error. The corresponding weak form of ˜ 1 (Ω) satisfying (1.6) reads: find ui ∈ W 2 ( A1 (u1 , v1 ) = (f1 , v1 ), ˜ 21 (Ω), ∀vi ∈ W (3.1) A2 (u2 , v2 ) = (f2 (x, u1 , Du1 ), v2 ), where A1 (u1 , v1 ) = A1 (u1 , v1 ) ≡ (a1 ∇u1 , ∇v1 ) + (b1 (x) · ∇u1 , v1 ) + (c1 u1 , v1 ), A2 (u2 , v2 ) = A2 (u2 , v2 ) ≡ (a2 ∇u2 , ∇v2 ) + (b2 (x) · ∇u2 , v2 ) + (c2 u2 , v2 ),
6
V. CAREY, D. ESTEP, AND S. TAVENER
˜ m (Ω) is the subspace of W m (Ω) are assumed to be coercive bilinear forms on Ω and W p p with zero trace on ∂Ω. We suppress the “cross” dependence on the other solutions except in a few remarks below. After introducing (conforming) discretizations Sh,i (Ω), we solve the discretized system ( A1 (U1 , χ1 ) = (f1 , χ1 ), ∀χi ∈ Sh,i (Ω). (3.2) A2 (U2 , χ2 ) = (f2 (x, U1 , DU1 ), χ2 ), In general, however, Sh,1 ( Sh,2 (or vice-versa) on Ω, and we may be forced to work with either Π1→2 f2 (U1 ) or more generally with f2 (x, Π1→2 U1 , Π1→2 DU1 ), where ˜ 1 (Ωi ) Πi→j is some projection from Sh,i to Sh,j . If the projection is to Sh,i from W 2 then we simply write the projection as Πi . The resulting discrete system becomes ( A1 (U1 , χ1 ) = (f1 , χ1 ), ∀χi ∈ Sh,i (Ω). (3.3) A2 (U2 , χ2 ) = (f2 (x, Π1→2 U1 , Π1→2 DU1 ), χ2 ), Primary Adjoint Problem We seek the error in a quantity of interest representable by a linear functional of the error e2 , where ui − Ui = ei denotes the pointwise errors. Note that a quantity of interest involving only u1 can be computed without solving for u2 , hence there is no loss of generality. The global adjoint problem, defined relative to the quantity of interest, is ( (1) (1) (1) (1) −∇ · a1 ∇φ1 − div(b1 φ1 ) + c1 φ1 + Lf2 (u1 )φ2 = 0, (1) (1) (1) (1) −∇ · a2 ∇φ2 − div(b2 φ2 ) + c2 φ2 = ψ2 , where Z
1
Lf2 (u1 )(u1 − U1 ) = 0 (1)
∂f2 (u1 s + U1 (1 − s)) ds ∂u1
(1)
is a linearization of f2 and φ1 and φ2 satisfy homogeneous Dirichlet boundary conditions. The corresponding weak formulation is ( (1) (1) A∗1 (φ1 , v1 ) + (Lf2 (u1 )φ2 , v1 ) = 0, ˜ 21 (Ω), ∀vi ∈ W (3.4) (1) (1) A∗2 (φ2 , v2 ) = (ψ2 , v2 ), where ( (1) (1) (1) (1) A∗1 (φ1 , v1 ) = (a1 ∇φ1 , ∇v1 ) − (div(b1 φ1 ), v1 ) + (c1 φ1 , v1 ) (1) (1) (1) (1) A∗2 (φ2 , v2 ) = (a2 ∇φ2 , ∇v2 ) − (div(b2 φ2 ), v2 ) + (c2 φ2 , v2 ).
(3.5)
Using the standard argument, we have the following error representation formula, ¡ (1) ¢ ¡ (1) ¢ ¡ (1) ¢ ¡ ¡ (1) ¢ (1) ¢ ψ , e = ψ2 , e2 = A∗2 φ2 , e2 = f2 (x, u1 , Du1 ), φ2 − A2 U2 , φ2 . (3.6) (1)
Observe that φ1 does not appear in the error representation formula. We define the primary adjoint problem as
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
(1)
7
(1)
˜ 21 (Ω). A∗2 (φ2 , v2 ) = (ψ2 , v2 ), ∀v2 ∈ W Remark 3.1. At first glance, it appears that we only need to solve the second adjoint equation and thus do not need to construct the linearization Lf2 . However, as seen in the linear algebra example, the analysis takes into account the transfer of error from the solution of the first component. Estimating this transferred error uses a nonlinear functional of the error to form the right hand sides in the “transfer adjoint problems”(2.3) and (3.11). We approximate this nonlinear functional using the linearization Lf2 . We evaluate the linearization at the computed solution U , which can be justified by using Taylor’s theorem and assuming that the error u − U is sufficiently small. (1) Adding and subtracting the projection of φ2 onto the primal approximation (1) space (Π2 φ2 ) in (3.6) yields (1)
(1)
(1)
(ψ2 , e2 ) = (f2 (x, u1 , Du1 ), (I − Π2 )φ2 ) − A2 (U2 , (I − Π2 )φ2 ) (1)
(1)
+ (f2 (x, u1 , Du1 ), Π2 φ2 ) − A2 (U2 , Π2 φ2 ). (3.7) To simplify later constructions, we introduce the notion of the weak residual of a solution component, namely Ri (Ui , χ; ν) = (fi (ν), χ) − Ai (Ui , χ; ν), and using this notation write (3.6) as (1)
(ψ (1) , e) = R2 (U2 , φ2 ; u1 ), indicating that this estimate depends on the solution u1 . 3.1. Transfer Error Analysis. The error representation (3.7) is not com¡ (1) ¢ putable since u1 is unknown. We add and subtract f2 (x, U1 , DU1 ), (I − Π2 )φ2 from the error representation formula (3.7) and use the definition of the approximate weak statement (3.2) to obtain ¡ (1) (1) ¢ (1) (ψ2 , e2 ) = f2 (x, U1 , DU1 ), (I − Π2 )φ2 − A2 (U2 , (I − Π2 )φ2 ) ¡ (1) ¢ + f2 (x, u1 , Du1 ) − f2 (x, U1 , DU1 ), φ2 ¢ ¡ ¡ (1) ¢ (1) = R2 U2 , (I − Π2 )φ2 ; U1 + f2 (x, u1 , Du1 ) − f2 (x, U1 , DU1 ), φ2 (3.8) The first term on the right of (3.8) is a traditional dual-weighted residual expression for the error arising from discretization of the second component while the remaining difference represents the transfer error that arises from using an approximation of u1 in defining the coefficients in the equation for u2 . The goal now is to estimate this transfer error and its effect on the quantity of interest. As with the linear algebra example in Section 2, we recognize the transfer error expression as a functional of error in u1 and define (1)
(f2 (x, u1 , Du1 ) − f2 (x, U1 , DU1 ), φ2 ) as a new quantity of interest. Then, we construct a secondary adjoint problem to compute the transfer error. In order to obtain a linear functional when f2 is nonlinear
8
V. CAREY, D. ESTEP, AND S. TAVENER
in u1 , we linearize, f2 (u1 ) ≈ f2 (U1 ) + Df2 (U1 ) × (u1 − U1 ), where Df is the Fr´echet derivative of f2 at U1 . The transfer error term becomes ¡
(1) ¢
Df2 (U1 ) × e1 , φ2
,
(3.9)
which is a linear functional of the error e1 that describes the effect of errors in U1 on the quantity of interest. Note that the Riesz representation theorem guarantees (2) (2) (2) the existence of a ψ1 such that (ψ1 , e1 ) equals (3.9), though ψ1 is not needed to evaluate the functional or compute the corresponding adjoint solution. Transfer Error Adjoint Problem To estimate the new quantity of interest, we define (¡ ¢ ¡ ¢ ¡ (2) ¢ (2) (2) (2) (2) a1 ∇φ1 , ∇v1 − div(b1 φ1 ), v1 + c1 φ1 , v1 ) + (Lf2 (u1 )φ2 , v1 = ψ1 , ¡ ¢ ¡ ¢ ¡ ¢ (2) (2) (2) a2 ∇φ2 , ∇v2 − div(b2 φ2 ), v2 + c2 φ2 , v2 = 0, ˜ 21 (Ω). ∀vi ∈ W The second equation has the trivial solution and the secondary adjoint problem reduces to the “transfer error adjoint problem” ¢ ¡ (2) ¢ ¡ (2) ¢ ¢ ¡ ¡ (2) (2) ˜ 21 (Ω). (3.10) ∀v1 ∈ W a1 ∇φ1 , ∇v1 − (div b1 φ1 , v1 + c1 φ1 , v1 = ψ1 , v1
The transfer error representation formula is given by ¡
¢ (2) (2) (2) ψ1 , e1 = A1∗ (φ1 , e1 ) = A1 (e1 , φ1 ) ¡ (2) (2) ¢ = f1 , (I − Π1 )φ1 − A1 (U1 , (I − Π1 )φ1 ),
(3.11)
where we have used Galerkin orthogonality to introduce the projection of φ onto the discretization space (as f1 does not depend on u). Inserting (3.11) into (3.8) yields ¡ ¢ ¡ (1) ¢ (1) ψ, e = f2 (x, U1 , DU1 ), (I − Π2 )φ2 − A2 (U2 , (I − Π2 )φ2 ) ¢ ¡ (2) (2) + f1 , (I − Π1 )φ1 − A1 (U1 , (I − Π1 )φ1 ). (3.12) or ¡ ¢ (1) (2) ψ, e = R2 (U2 , (I − Π2 )φ2 ; U1 ) + R1 (U1 , (I − Π1 )φ1 ). Remark 3.2. If the model problem includes coupling in the coefficients of the second differential operator, i.e., x ∈ Ω, −∇ · a1 (x)∇u1 + b1 (x) · ∇u1 + c1 (x)u1 = f1 (x), −∇ · a2 (x, u1 )∇u2 + b2 (x, u1 ) · ∇u2 + c2 (x, u1 )u2 = f2 (x, u1 , Du1 ), x ∈ Ω u1 = u2 = 0, x ∈ ∂Ω, (3.13) then the error representation formula for a quantity of interest that depends on u2 alone is ¡ ¢ (1) ψ, e = R2 (U2 , (I − Π2 )φ2 ; u1 ).
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
9
Since this is not computable, we replace each term in the weak residual with the same term evaluated at U1 ; yielding ¡
¢ ¡ (1) (1) ¢ ψ, e = R2 (U2 , (I − Π2 )φ2 ; U1 ) + f2 (u1 ) − f2 (U1 ), φ2 ¡ ¢ ¡ ¢ − (a2 (u1 ) − a2 (U1 ))U2 , φ(1) − (b2 (u1 ) − b2 (U1 )) · ∇U2 , φ(1) ¡ ¢ − (c2 (u1 ) − c2 (U1 ))U2 , φ(1)
We linearize f2 , a2 , b2 , and c2 around U1 to obtain an approximate transfer error term ¡ ¡ ¡ ¡ (1) ¢ (1) ¢ (1) ¢ (1) ¢ Df2 (e1 ), φ2 + Da2 (e1 )∇U2 , ∇φ2 + Db2 (e1 ) · ∇U2 , φ2 + Dc2 (e1 )U2 , φ2 . This is a linear functional on L2 (Ω), which we use as data to define the “transfer” error adjoint problem and derive a corresponding a posteriori error representation. For details on how to compute a quantity of interest that depends on u1 and u2 (so that the choice of linearizations for the coefficients in the equation for u2 enter directly into the “primary” error contribution) see [9]. Remark 3.3. For a “lower triangular” one-way coupled system of N elliptic equations and a quantity of interest based on the N th component, we solve N total “single physics” adjoint problems and construct the error representation ¡
N ¢ X ¡ ¢ ψN , eN = RN −i+1 UN −i+1 , φ(i) ; U . i=1
We then solve a sequence of adjoint problems, as the corresponding linear functional for the ith adjoint problem (i > 1) can be defined recursively (assuming the coupling only occurs through the right-hand side) as ¯ ¶ i−1 µ X DfN +1−j ¯¯ (j) (e ), φ . i Dui ¯ j=1
U
This extends to coupling in all the coefficients as above. 3.2. Numerical examples. The following three numerical examples highlight the features of the analysis and the importance of accounting for the transfer error. In the following computations, we approximately solve all adjoint problems using continuous, piecewise quadratic elements in order to be able to evaluate the interpolants arising from Galerkin orthogonality. We denote these approximate adjoints solutions by Φ and use them in place of φ in the error representation (3.12). For adaptive mesh refinement, we write the estimate as a sum of element contributions and derive a bound by introducing norms. We base the adaptive mesh refinement on the standard optimization approach using the Principle of Equidistribution [6] applied to the bound. We refine elements whose element contribution to the error bound is greater than half a standard deviation from the mean error contribution or refine a fixed fraction of the elements with the greatest element contributions, whichever criterion yields the greater refinement. We do not do any mesh coarsening, smoothing, or edge flips. 3.2.1. Example. This example demonstrates the fact that the transfer error can be significant even if the individual components u1 and u2 are well resolved. We
10
V. CAREY, D. ESTEP, AND S. TAVENER
consider a simple system −∆u1 = sin(4πx) sin(πy), −∆u2 = b · ∇u1 , u = 0, where 2 b= π
µ
25 sin(4πx) sin(πx)
¶
µ ,
f (u) =
(x, y) ∈ Ω, (x, y) ∈ Ω, (x, y) ∈ ∂Ω,
sin(4πx) sin(πy) b · ∇u1
(3.14)
¶ Ω = ([0, 1], [0, 1]).
The quantity of interest is the solution value of u2 at (.25, .25), which we estimate using a smooth delta function approximation with localized support. The corresponding global adjoint problem is (1) (1) −∆φ1 + Lf (u1 )φ2 = 0, (x, y) ∈ Ω (1) (3.15) −∆φ2 = δxreg (x, y) ∈ Ω, ˜ , φ = 0, (x, y) ∈ ∂Ω, where δxreg is a regularized delta function and x ˜ = (.25, .25). Our primary adjoint ˜ problem is (1)
−∆φ2 = δxreg ˜ ,
(x, y) ∈ Ω,
φ = 0, (x, y) ∈ ∂Ω,
The secondary adjoint problem is ( (1) (2) ∆φ1 = ∇ · (bφ2 ), (x, y) ∈ Ω, φ(2) = 0, (x, y) ∈ ∂Ω.
(3.16)
The primal system was solved using identical standard continuous piecewise linear finite element discretizations for u1 and u2 . We plot the results in Fig. 3.1 and show (2) the error contributions in Table 3.1. While the adjoint solution Φ1 in Fig. 3.1(d) is concentrated near the location of the quantity of interest, it has non-trivial spatial structure and the transfer error represents 14% of the total error. Primary Error 0.0042
Transfer Error 0.0006
Table 3.1: Error contributions for Example 3.2.1
3.2.2. Example. This example illustrates the importance of computing the transfer error, since for this problem simply forcing the “primary” error contribution to be small (by refining the second mesh only) does not provide any accuracy in the desired quantity of interest. We reconsider (3.14) but with quantity of interest equal to the average value of u2 over the whole domain. The exact solution has zero average value on Ω. We solve both components of the primary problem on an identical coarse initial mesh, but only adapt and refine the mesh for u2 using the traditional weighted residual, the first “primary” error term in (3.12) while neglecting the second “transfer error” term in (3.12). We show the results in Fig. 3.2 and Table 3.2.
11
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
1 0.5
0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 1
0 −0.5 −1 −1.5 1 0.5 0 0
0.4
0.2
1
0.8
0.6
0.5 0.2
0 0
(a) U1
0.4
0.6
1
0.8
(b) U2
2 0.2 0.15 0.1 0.05 0 −0.05 −0.1 1
1.5 1 0.5 0 1 0.5 0.2
0 0
0.4
0.6
1
0.8
0.5 0.2
0 0
(1)
0.4
0.6
1
0.8
(2)
(c) Φ2
(d) Φ1
Fig. 3.1: Example 3.2.1: Primary and (non-zero) adjoint solutions computed on uniformly fine meshes. The adjoint solutions are largest near the region of the quantity of interest u2 (.25, .25)
0.4 0.2 0.5 0
0
−0.2
−0.5
−0.4
−1
−0.6 1
−1.5 −2 1 0.5 0.5 0
0
0.2
0.4
(a) U1
0.6
0.8
1
0 0
0.2
0.4
0.6
0.8
1
(b) U2
Fig. 3.2: Example 3.2.2: Adaptivity based on the standard discretization error estimate for the “primary” error, ignoring the “transfer” error. Only the mesh for U2 is refined.
12
V. CAREY, D. ESTEP, AND S. TAVENER
Primary Error 0.00005
Transfer Error 0.110
Table 3.2: Error contributions for Example 3.2.2
Ignoring the transfer error and the implied need to refine the first component produces a completely unsuitable adaptive procedure. It is clear from Fig. 3.2 that the average value of the second component is far from zero, and the actual computational value is −0.2245. The estimated transfer error of 0.1 is in fact an underestimate since (2) Φ1 is based on the highly inaccurate solution U1 which is computed on a very coarse mesh. The transfer error dominates the computation and this error cannot be reduced without refining the mesh for u1 . 3.2.3. Example. The third example shows that an “optimal” adaptive mesh for the quantity of interest that depends only on u2 may actually involve a richer discretization of u1 than u2 . We consider the system (3.14) with the quantity of interest equal to the average value of u2 over the whole domain and initial coarse meshes as in the previous example, but use the transfer error contribution to adapt the mesh for u1 and the primary contribution to adapt the mesh for u2 so that the total error is less than 10−4 . The resulting meshes are shown in Fig. 3.3, and illustrate that despite the fact that the quantity of interest only involves u2 , the error inherited from u1 is the most important contribution to consider. In this problem, the strong influence of the transfer error is a result of the dependence of u2 on the gradient of u1 , which a priori has lower order accuracy. Similar behavior could also arise when u2 just depends on u1 .
0.8 0.6 1
0.4
0.5
0.2 0
0 −0.2 −0.5
−0.4 −0.6
−1 1 0.5 0 0
(a) U1
0.2
0.4
0.6
0.8
1
−0.8 1 0.5 0 0
0.2
0.4
0.6
0.8
1
(b) U2
Fig. 3.3: Example 3.2.3: Adaptivity based on the full estimate that accounts for “primary” and “transfer” errors. The quantity of interest U2 (.25, .25) is more sensitive to errors in U1 than U2 .
4. Interpolation error analysis. We use a multiscale discretization for the “fully” adaptive example 3.2.3, i.e., the components u1 and u2 were computed on different meshes, see Figure 3.3. This raises the issue of understanding the effect of translating one component onto the mesh of the other component when performing the
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
13
integration necessary to form the discrete equations. Integration involving functions defined on different meshes these quantities may be complicated, as illustrated in Figure 4.1. Mesh for U1
Mesh for U2
Fig. 4.1: The problem of translation between meshes. Finite element functions on one mesh are generally not smooth on another mesh.
In particular, traditional quadrature formulae based on sets of specific points may not preserve the accuracy required for effective computation because a function defined on a different mesh is generally not sufficiently smooth. For example, the integrand (f2 , χ) is piecewise discontinuous on every element τi of mesh 2 in Example 3.2.3 as b·∇U1 is continuous only within elements of the mesh for U1 . In general, if the meshes are not congruent, the integrand is C 0 at best. Using a “traditional” higher order quadrature rule will not necessarily lead to the expected increase in accuracy as the integrand (f2 , χ) does not have sufficient regularity. Possible solutions include either the determination of local intersections of simplices and/or hexahedra or the construction of a global union mesh. However, both solutions are computationally expensive and the global solution often requires several times more memory than the storage of the two individual meshes, especially for three dimensional problems. 4.1. Projections from mesh 1 to mesh 2. Instead of constructing a union mesh, we use a projection Π1→2 from S1,h to S2,h and solve the discrete system given by (3.3). This introduces additional sources of error. Starting from the error representation formula (3.6), we add and subtract (1)
f2 (x, Π1→2 U1 , Π1→2 DU1 , (I − Π2 )φ2 ), yielding (1)
(ψ (1) , e) = (ψ2 , e2 ) ¡ ¡ (1) ¢ (1) ¢ = f2 (x, Π1→2 U1 , Π1→2 DU1 ), (I − Π2 )φ2 − A2 U2 , (I − Π2 )φ2 ¡ (1) ¢ + f2 (x, u1 , Du1 ) − f2 (x, Π1→2 U1 , Π1→2 DU1 ), φ2 . ¡ (1) ¢ Adding and subtracting f2 (x, U1 , DU1 ), φ2 produces ¡ ¡ (1) ¢ (1) ¢ (ψ (1) , e) = f2 (x, Π1→2 U1 , Π1→2 DU1 ), (I − Π2 )φ2 − A2 U2 , (I − Π2 )φ2 ¡ (1) ¢ + f2 (x, u1 , Du1 ) − f2 (x, U1 , DU1 ), φ2 ¡ (1) ¢ + f2 (x, U1 , DU1 ) − f2 (x, Π1→2 U1 , Π1→2 DU1 ), φ2 .
14
V. CAREY, D. ESTEP, AND S. TAVENER
The first two terms on the right represent the primary discretization error for a functional of the the second component, the third term on the right represents the transfer error (3.11), and the fourth term is a new expression that represents the error from the projection Π1→2 . The projection error can be decomposed as ¡ (1) ¢ Π1→2 f2 (x, U1 , DU1 ) − f2 (x, Π1→2 U1 , Π1→2 DU1 ), φ2 ¡ (1) ¢ + (I − Π1→2 ) f2 (x, U1 , DU1 ), φ2 (4.1) The first inner product in (4.1) can be computed (with some effort) on Ω2,h . However, computing the second term raises the same numerical issues that caused the adoption of the projection Π1→2 in the first place! We handle this term using the Monte Carlo techniques described in Section 4.3. 4.2. Projections from mesh 2 to mesh 1. Complications from the use of projections also arise in computations with the solution of the secondary adjoint problem. (1) The secondary adjoint problem domain is Ω1,h but φ2 is computed naturally on Ω2,h . The new error representation formula for the transfer error becomes ¡ ¡ (1) ¢ (1) ¢ Df2 (U1 ) × e1 , Π2→1 φ2 + Df2 (U1 ) × e1 , (I − Π2→1 )φ2 which is the error contribution arising from the transfer as well as an additional term that is large when the approximation spaces are significantly different. For example, this term is important when the original system is multiscale. The implicit ψ (2) for the transfer error adjoint is now ¡ ¡ ¡ (2) ¢ (1) ¢ (1) ¢ f2 (u1 ) − f2 (U1 ), φ2 = Df2 (U1 ) × e1 , Π2→1 φ2 = ψ1 , e1 ¡ (1) ¢ The additional term Df2 (U1 ) × e1 , (I − Π2→1 )φ2 is a linear functional, so we may define an additional “tertiary” adjoint problem to estimate this quantity. Projection{“Tertiary”} Error Adjoint Problem This problem has the same form as the transfer error adjoint (3.10), but with (3) data ψ1 satisfying ¡ (3) ¢ ¡ (1) ¢ ψ1 , e1 = Df2 (U1 ) × e1 , (I − Π2→1 )φ2 . The resulting error representation formula is ¡ (3) ¢ ¡ ¡ (3) (3) ¢ (3) ¢ ψ1 , e1 = f1 , (I − Π1 )φ1 ) − A1 (U1 , (I − Π1 )φ1 = R1 , (I − Π2→1 )φ1 (4.2) The error representation is therefore (1)
(1)
(2)
(3)
(ψ2 , e2 ) = R2 (U2 , (I − Π2 )φ2 ; U1 ) + R1 (U1 , (I − Π1 )(φ1 + φ1 )) (4.3) ¡ ¡ (1) ¢ (1) ¢ + Π1→2 f2 (U1 ) − f2 (Π1→2 U1 ), φ2 + (I − Π1→2 )f2 (U1 ), φ2 . Remark 4.1. Traditional simplex-based numerical integration methods that interrogates U1 at cubature points can be thought of as projecting the integrand f (U1 )χ2 into a specific polynomial space Pτ defined on each simplex τ of the mesh for U2 and then integrating exactly. We may express this “cubature error” as a projection error and construct a corresponding error representation formula in a similar manner. Cubature error resulting from the fact that integration was not performed on a “union” mesh of two piecewise polynomial spaces may always be viewed as projection error.
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
15
Remark 4.2. In this discussion, we assume that the adjoint problems are solved using approximations spaces that are compatible with the corresponding primal approximation space, e.g., using higher-order Lagrange elements on the same mesh. In practice, different meshes may be used for the primal and adjoint solves. However, this introduces new projection operators between the corresponding approximation spaces as well as the additional terms due to the loss of Galerkin orthogonality. We confine ourselves to merely alluding to the notational complexities and length of the resulting error representation. 4.3. Monte Carlo Integration. Interpolation-based projections suffer from mesh-aliasing difficulties. An extreme example is given in Fig. 4.2. 1
U1 Π1→2U1
0
−1
0
0.2
0.4
0.6
0.8
1
Fig. 4.2: Interpolation errors for two meshes on Ω = [0, 1].
For a more practical example, we construct two quasi-uniform, unstructured meshes, 1 and 2, both of size h on Ω = [0, 1] × [0, 1], and take the piecewise linear interpolantR fI1 of the function R f = sin(20hx) sin(20hy) on mesh 1. We first compute Ie = Ω f dx and I1 = Ω fI1 dx exactly and then construct three different approximations Igauss , IΠ and ISamp as follows. 1. Igauss : Using a 3rd order, 4-point quadrature rule [16] on the triangles of mesh 2 by interpolating fI1 at the corresponding quadrature points. 2. IΠ : Projecting fI1 on to mesh 2 by interpolating fI1 at the nodes of mesh 2 and then using exact integration. 3. ISamp : Performing the integration via a uniform weight quadrature rule using the quadrature points corresponding to the 4-point quadrature rule employed by Igauss . We show the accuracy in Table 4.1. Note that the work for all 3 methods is roughly the same. The smallest of the projection errors |I1 −Igauss |, is larger than the interpolation |Ie − I1 | 0.000187
|I1 − Igauss | 0.000246
|I1 − IΠ | 0.0060
|I1 − ISamp | 0.00041
Table 4.1: Errors in various approximations of Ie . error |Ie − I1 |. The error in |I1 − IΠ | is a factor of 10 larger than |I1 − Igauss | and |I1 − ISamp |, which, for this problem, amounts to a factor of h−1 . Note that the 4-point gauss quadrature rule is only slightly more accurate than the sampling rule ISamp . Motivated by the example, we employ pseudo-random Monte Carlo integration using p random uniformly-distributed sample points on the reference element. The
16
V. CAREY, D. ESTEP, AND S. TAVENER
main difficulty (and computational expense) when integrating on Ω2,h is the evaluation of U1 at each random sample point since this involves locating the point in the appropriate element in Ω1,h . Nominally, this process requires (O(N )) operations per sample point, where N is the number of degrees of freedom for U1 , hence O(M N ) operations for the integration, where M is the number of degrees of freedom for U2 . However, this approach may be greatly accelerated by using a geometric implementation of the assembly and point search algorithms.
τi ∈ Ω2,h pi1 pi2 pi+1 1 u1,h
Fig. 4.3: Monte Carlo integration point search
We illustrate the search algorithm in Fig. 4.3. We generate a random integration point p11 in τ1 ∈ Ω2,h , and determine the containing element of Ω1,h . This could potentially involve a full search of Ω1,h but as this is the initial element, a good starting guess for element location could be provided as an input. Once a matching simplex is found in Ω1,h , the computation is performed and the next integration point p12 is generated. Moreover, the last matching simplex is stored, so the geometric search using edge/face neighbors and barycentric coordinates to guide neighbor selection for the next point is very fast. When the integration is finished, we select the next element to be an edge/face neighbor. Now when we generate p21 , we have a good starting point, namely the last match in p1S which should be “close” to the real element containing p21 . The assembly routine keeps selecting edge neighbors until it has looped over all elements recursively. This algorithm works even with a primitive data structure as long as recursion is employed. If the number of mesh elements is large, however, this may not be practical due to recursion limits. A non-recursive algorithm could lead to termination before all element contributions for the mesh were calculated, as the next element returned by the search could have all edge/face neighbors whose element contributions had already been calculated. The algorithm would have to “restart” from an element that has not been computed. On quasi-uniform meshes with no fine scale features in the geometry, the number of “restarts” also grows logarithmically with the number of elements. Of course, with a more sophisticated data structure, either octree based or, for example, a mesh where the elements had been ordered by the use of a space filling curve, the need for restarting would be eliminated. When the meshes for Ω1,h and Ω2,h are both quasi-uniform on Ω, the number of elements tested in Ω is bounded by some h-independent constant for each integration point. Obviously this is not the case for general adapted or anisotropic meshes, but in practice the number of searches grows at most logarithmically with the numbers of degrees of freedom in u1 . The convergence of this Monte Carlo integration scheme follows from standard results (see [11]) as the integrand can always be defined as the
17
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
sum of integrals of continuous functions on individual simplices of the union mesh of Ω1,h and Ω2,h . 4.4. Numerical examples. We demonstrate the significance of the projection errors with two examples. 4.4.1. Example. The first example illustrates how the projection error can influence a typical computation. We consider a system defined by (3.14), with two randomly-generated initial meshes for u1 and u2 . The initial mesh for u1 is finer than for u2 in order to reduce the transfer error. The quantity of interest in this computation is the average value of u2 . We show the results in Fig. 4.4 and Table 4.2. We use a local projector Π1→2,τ given by interpolation at the Gauss points (third order 3 point simplex rule) of simplices τ in Sh,2 . Use of this projector would integrate (U1 , U2 ) exactly if the meshes were identical. The solution using this projector is given by Fig. 4.4(b). This is compared against a 16 point Monte Carlo computation illustrated by Fig. 4.4(c).
2 1.5
1
1 0.5
0.5 0
0
−0.5 −0.5
1 0.8
−1 1
0.6 1
0.4
0.8
0.5
0.6
0.2
0.4 0
0.2
0
0
(a) U1
0.2
0
0.4
0.6
0.8
1
(b) U2 computed with Π1→2,τ
0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 1
0.5
0
0
0.2
0.4
0.6
0.8
1
(c) U2 computed with high-sample Monte Carlo integration
Fig. 4.4: Example 4.4.1: The role of projection errors in non-aligned meshes. Note that the magnitude and oscillation of U2 computed with Π1→2,τ , shown in (b) are incorrect (a fine scale U2 is given by Fig. 3.1).
18
V. CAREY, D. ESTEP, AND S. TAVENER
Primary Error 0.003533
Transfer Error 0.021589
Projection Error 0.007908
Table 4.2: Example 4.4.1: Error contributions for computation shown in Fig. 4.4.
2 2 As discussed in Section 4.2, projection from Sh,2 to Sh,1 can also lead to significant inaccuracies in computing the transfer error, necessitating the computation of a tertiary adjoint problem (4.2). 2 2 to Sh,1 can 4.4.2. Example. As discussed in Section 4.2, projection from Sh,2 also lead to significant inaccuracies in computing the transfer error, necessitating the computation of a tertiary adjoint problem (4.2). This example shows that computations with significant differences in mesh scale can contribute significantly to the error. We again use the system in Example 3.2.1 with quantity of interest the point value at (.15, 15), starting with a coarse identical initial mesh for u1 and u2 , but refining only the mesh for u2 . There is no projection error as Sh,2 ⊆ Sh,1 . However, when we compute the transfer error, we ignore the fact that a natural choice of decomposition for the computation is integration over the simplices of Sh,2 . Instead, we use (1) interpolation of φ2 at the quadrature points at the simplices of Sh,1 . To compute (I − Π), we employ the actual nesting of the two meshes to perform an accurate (up (3) to quadrature error on the fine scale mesh) computation of φ1 . We show the results in Table 4.3 and Fig. 4.5.
Primary Error 0.000713
Transfer Error 0.0905
Projection Error 0
Tertiary Error 0.0325
Table 4.3: Example 4.4.2: Error contributions for the computation shown in Fig. 4.5.
1 0.5 0 −0.5 1 −1 0
0.5 0.2
0.4
0.6
0.8
1 0
(3)
Fig. 4.5: Example 4.4.2: Tertiary adjoint solution Φ1 which estimates the projection error in computing the transfer error
5. An adaptive algorithm for the operator decomposition - finite element method. An adaptive algorithm that takes into account all the possible
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
19
sources of error is given below. while (the total error is less than TOL) do Compute U1 using standard integration. Compute U2 using 16 point M.C. integration for the coupling term. (1) Compute Φ2 using standard integration (2) Compute Φ1 for given adjoint data using 16 point M.C. integration if (the sum of two error contributions is greater than TOL) then Refine both meshes based on the primary error contributions for U2 and the transfer error contributions for U1 . else Compute the projection error by comparing with a 64 point M.C. integration. (3) Compute Φ1 . if (the total error is greater than TOL) then Refine both meshes based on the primary and projection error contributions for U2 and the transfer and tertiary error contributions for U1 . end if end if end while The algorithm drives the primary and transfer error contributions to within a specified error tolerance, and then checks for projection error by using 64 point Monte Carlo integration as an approximation to the identity operator I in (4.3) and attempts to correct the projection error by refinement as well. Any projector could be substi(2) tuted for the M.C. integration used in computing U2 and Φ1 . We select the use of 16 sample points for the Monte Carlo integration based on our experience from a series of numerical experiments where different functions were interpolated on a quasi uniform mesh, and integrated. This interpolant was then integrated using Monte Carlo with 2N sample points per simplex on a different quasiuniform mesh (with the same approximate h); N = 4 gave the best tradeoff between speed and accuracy. 5.1. Examples. We describe two applications of the algorithm to one way coupled systems using different meshes for each solution component. In both examples, we start with identical coarse initial meshes (quasi-uniform with h ≈ .125) and adapt each mesh until both the primary and transfer error formulas are less than 10−4 . We control projection error using Monte Carlo integration. 5.1.1. Example. In the first example, we approximate the value of u2 at (.25, .25), where (u1 , u2 ) solves 2 −∆u1 = 64π sin 4π(x − .75 + |x − .75|) sin 4π(y − .75 + |y − .75|), −∆u2 = u1 , u1 = u2 = 0, with Ω = ([0, 1], [0, 1]). The corresponding adjoint problem is ( (1) −∆φ2 = δreg (x0 ), (x, y) ∈ Ω, (1) φ2 = 0, (x, y) ∈ ∂Ω,
(x, y) ∈ Ω, (x, y) ∈ Ω, (x, y) ∈ ∂Ω,
20
V. CAREY, D. ESTEP, AND S. TAVENER
with x0 = (.25, .25). The transfer error adjoint problem is ( (2) (1) −∆φ1 = φ2 , (2) φ1 = 0,
(x, y) ∈ Ω, (x, y) ∈ ∂Ω.
1
1
0.5 0.5 0 1 0.8 0
0.6 0.4 0.2
−0.5 1
0.8
0.6
0.4
0.2
0 1
0.5
0
0
0
0.2
0.4
0.6
0.8
1
(a) U1 computed using a uniformly fine mesh, (b) U2 for the quantity of interest u2 (0.25, 0.25) axes rotated
0.8 0.6 0.4 0.2
1 0.8
0
0.6
−0.2
0.4
−0.4
0.2
−0.6
0
0
−0.4 1
1
0.2
−0.2 0.4
0.5
0.6 0.8
0.6
0.4
0.8 0.2
0
1
0 0
0.2
0.4
0.6
0.8
1
(c) U1 for the quantity of interest u2 (0.25, 0.25), (d) U1 for the quantity of interest u2 (0.9, 0.9) axes rotated
Fig. 5.1: Example 5.1.1: Example of computational efficiency: U1 may be computed on a coarse discretization yet U2 may be determined with sufficient accuracy.
The accurate computation of the quantity of interest u2 (0.25, 0.25) does not require the fine scale features of u1 near (0.75, 0.75) to be resolved. However, a quantity of interest equal to the value of u2 at (0.9, 0.9), near the localized features of u1 requires better resolution of the details of u1 . The adapted solutions U1 for both quantities of interest and given in Fig. 5.1(c) and Fig. 5.1(d) respectively. 5.1.2. Example. We now consider an example where convection in component u1 creates the need for refinement in u1 remote from the the goal-oriented refinement
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
in u2 .
3 −100kx−x0 k2 , x ∈ Ω, −∆u1 − b · ∇u1 = 10 e 3 −100kx−x1 k2 −∆u2 = 10 e u1 , x ∈ Ω, u1 = u2 = 0, x ∈ ∂Ω
21
(5.1)
where b = (100 40)> , x0 = (.75, .75), x1 = (.1, .5), and the quantity of interest is the point value u2 (x2 ), x2 = (.2, .5). The corresponding adjoint problem for the primary error contribution is ( (1) −∆φ2 = δreg (x2 ), x ∈ Ω, (1) φ2 = 0, x ∈ ∂Ω, while the corresponding transfer error adjoint problem is ( 2 (1) (2) (2) −∆φ1 + b · ∇φ1 = 103 e−100kx−x1 k φ2 , x ∈ Ω, (2) φ1 = 0, x ∈ ∂Ω. (2)
The adjoint solution Φ1 in Fig. 5.2(c) shows the influence of the convection term in the equation for u1 . When the quantity of interest is a value of u2 in the convective region of influence of the localized source term in the equation for u1 , the solution for u1 is resolved “upstream” of the location of the quantity of interest as shown in Fig. 5.2(a). When the quantity of interest is a value of u2 away from the convective region of influence of the localized source term in the equation for u1 , the adjoint solution (2) φ1 has an similar structure to that shown in Fig. 5.2(c) but has much smaller magnitude. The resulting mesh for U1 need not even be detailed enough to eliminate the numerical oscillation (from not satisfying the corresponding P´eclet mesh condition). This situation is illustrated by Fig. 5.2(d), where the choice of quantity of interest is u2 at (0.15, 0.15). 6. Conclusion. In this paper, we perform an a posteriori error analysis of a multiscale operator decomposition finite element method for the solution of a system of one-way coupled elliptic problems. The analysis specifically accounts for the effects arising from multiscale operator decomposition, including the issues: (1) errors in the solution of each component propagate into the solutions of the other components; and (2) transferring information between different representations potentially introduces new error. We estimate the various sources of errors by defining auxiliary adjoint problems whose data are related to errors in the information passed between components. Through a series of examples, we demonstrate the importance of accounting for the contributions to the error arising from multiscale operator decomposition. We also devise an adaptive discretization strategy based on the error estimates that specifically controls the effects arising from operator decomposition. Finally, we demonstrate the usefulness of Monte Carlo integration methods for dealing with mismatch between discretizations of different components. We extend this analysis to a “fully-coupled” system in the form of (1.4) in Part II of this paper [3]. We address the important issue that the adjoint operator associated with the fully coupled system and an operator decomposition solution are not generally equal. This difference requires additional strategies for error control. We consider the use of non-interpolatory projectors based on averaging to reduce both transfer and projection error in [4].
22
V. CAREY, D. ESTEP, AND S. TAVENER
1.5 5 1 4 0.5 3 0 2 −0.5 1
1 0 1 0.5 0 0
0.2
0.4
0.8
0.6
1
0.5 0 0
(a) U1 for quantity of interest u2 (.2, .5)
0.2
0.4
0.6
0.8
1
(b) U2 for quantity of interest u2 (.2, .5)
0.5 1.5
0
1
−0.5 1
0.5
0.8
0
0.6 −0.5 1
0.4
1
0.2
0.8
0.5
0.6 0.4
0
0 (2)
(c) Φ1
0.2
0.4
0.6
0.8
1
for quantity of interest u2 (.2, .5)
0
0.2 0
(d) U1 for quantity of interest u2 (.15, .15)
Fig. 5.2: The role of convection in Example 5.1.2. Note that altering the location of the quantity of interest alters the density and location of the resulting adapted meshes (same adaptive criteria).
REFERENCES [1] Wolfgang Bangerth and Rolf Rannacher, Adaptive Finite Element Methods for Differential Equations, Birkhauser Verlag, 2003. [2] R. Becker and R. Rannacher, An optimal control approach to a posteriori error estimation in finite element methods, Acta Numerica, (2001), pp. 1–102. [3] Varis Carey, Donald Estep, and Simon Tavener, A Posteriori analysis and adaptive error control for operator decomposition methods for elliptic systems II: Fully coupled systems, In preparation, (2007). [4] , Averaging based projections in operator decomposition methods for elliptic systems, In preparation, (2007). [5] K. Eriksson, D. Estep, P. Hansbo, and C. Johnson, Introduction to adaptive methods for differential equations, in Acta numerica, 1995, Acta Numer., Cambridge Univ. Press, Cambridge, 1995, pp. 105–158. [6] , Computational Differential Equations, Cambridge University Press, Cambridge, 1996. [7] D. Estep, A posteriori error bounds and global error control for approximation of ordinary differential equations, SIAM J. Numer. Anal., 32 (1995), pp. 1–48. [8] D. Estep, M. Holst, and M. Larson, Generalized Green’s functions and the effective domain of influence, SIAM J. Sci. Comput., 26 (2005), pp. 1314–1339. [9] D. Estep, M. G. Larson, and R. D. Williams, Estimating the error of numerical solutions of
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
23
systems of reaction-diffusion equations, Mem. Amer. Math. Soc., 146 (2000), pp. viii+109. [10] D. Estep, S. Tavener, and T. Wildey, A posteriori analysis and improved accuracy for an operator decomposition solution of a conjugate heat transfer problem. SINUM, in revision, 2006. [11] George S. Fishman, Monte Carlo, Springer Series in Operations Research, Springer-Verlag, New York, 1996. Concepts, algorithms, and applications. ¨ li, Adjoint methods for PDEs: A posteriori error analysis and postpro[12] M. Giles and E. Su cessing by duality, Acta Numerica, (2002), pp. 145–236. [13] V. Ginting, D. Estep, J. Shadid, and S. Tavener, An a posteriori analysis of operator splitting. to appear in SINUM, 2006. [14] G.I. Marchuk, On the theory of the splitting-up method, in Proceedings of the Second Symposium on Numerical Solution of Partial Differential Equations, SVNSPADE, Academic Press, New York, 1970, pp. 469–500. [15] , Splitting and alternating direction methods, in Handbook of Numerical Analysis, vol. I, North-Holland, New York, 1990, pp. 197–462. [16] Gilbert Strang and George J. Fix, An Analysis of the Finite Element Method, Prentice-Hall Inc., Englewood Cliffs, N. J., 1973. Prentice-Hall Series in Automatic Computation.
1. Introduction. Multiscale operator decomposition is a widely used technique for solving multiphysics, multiscale problems. The general approach is to decompose the multiphysics problem into components involving simpler physics over a relatively limited range of scales, and then to seek the solution of the entire system through some sort of iterative procedure involving solutions of the individual components. This approach is appealing because there is generally a good understanding of how to solve a broad spectrum of single physics problems accurately and efficiently, and because it provides an alternative to accommodating multiple scales in one discretization. However, multiscale operator decomposition presents an entirely new set of accuracy and stability issues, some of which are obvious and some subtle, and all of which are difficult to correct. We motivate multiscale operator decomposition for elliptic systems by considering a model of a thermal actuator. A thermal actuator is a MEMS (microelectronic mechanical switch) device. A contact rests on thin braces composed of a conducting material. When a current is passed through the braces, they heat up and consequently expand to close the contact. The system is modelled by a system of three coupled equations, each representing a distinct physical process. They are an electrostatic ∗ Department of Mathematics, Colorado State University, Fort Collins, CO 80523 ([email protected]). V. Carey’s work is supported in part by the Department of Energy (DE-FG02-04ER25620) † Department of Mathematics and Department of Statistics, Colorado State University, Fort Collins, CO 80523 ([email protected]). D. Estep’s work is supported in part by the Department of Energy (DE-FG02-04ER25620, DE-FG02-05ER25699, DE-FC02-07ER54909), the National Aeronautics and Space Administration (NNG04GH63G), the National Science Foundation (DMS0107832, DMS-0715135, DGE-0221595003, MSPA-CSE-0434354, ECCS-0700559), Idaho National Laboratory (00069249), and the Sandia Corporation (PO299784) ‡ Department of Mathematics, Colorado State University, Fort Collins, CO 80523 ([email protected]). S. Tavener’s work is supported in part by the Department of Energy (DE-FG02-04ER25620)
1
2
V. CAREY, D. ESTEP, AND S. TAVENER Vsig
Vswitch
Fig. 1.1: Sketch of a thermal actuator.
current equation ∇ · (σ∇V ) = 0,
(1.1)
governing potential V (where current J = −σ∇V ), a steady-state energy equation ∇ · (κ(T )∇T ) = σ(∇V · ∇V ),
(1.2)
governing temperature T , and a linear elasticity equation giving the steady-state displacement d, ¡ ¢ ¡ ¢ ∇ · λ tr(E)I + 2µE − β(T − Tref )I = 0, E = ∇d + ∇d> /2. (1.3) Using multiscale operator decomposition, the complete system (1.1-1.3) is decomposed into three components, each of which is solved with a code specialized to the particular type of physics. Notice that the electric potential V can be calculated independently of T and d. The temperature T can be calculated once the electric potential V is known, while the calculation of displacement d requires prior knowledge of T , and therefore of V . In general, we can write a coupled elliptic system on a domain Ω in the form L (x, u1 , Du1 , . . . , un , Dun ) = 0 1 .. , x ∈ Ω. (1.4) . L (x, u , Du , . . . , u , Du ) = 0 n 1 1 n n A natural form of operator decomposition is to split the global multi-physics problem into n “single-physics” components that are solved individually. In general, the solution of each component requires knowledge of the solutions of all the other components; the full problem requires some form of iteration to obtain the solution. It is possible to impose conditions on the system, the components, and the coupling that allow for an a priori convergence analysis. However, operator decomposition is problematic in practice because it is very difficult to verify such conditions
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
3
and often impractical to satisfy them. Indeed, numerical solutions obtained via operator decomposition are affected significantly by the specific choice of decomposition. In this paper, we perform an a posteriori error analysis of a multiscale operator decomposition finite element method for the solution of a system of coupled elliptic problems. The components of the problem are solved in sequence using independent discretizations. The goal is to compute accurate computational error estimates that specifically account for the effects arising from operator decomposition. We also devise an adaptive discretization strategy based on the error estimates that controls the effects arising from multiscale operator decomposition. The a posteriori analysis in this paper is based on a well known approach involving variational analysis, residuals and the generalized Green’s function solving an adjoint problem [1, 2, 5–7, 9, 12]. However, we modify this approach to accommodate several new features arising from the operator decomposition. Three important issues addressed here are: (1) errors in the solution of each component propagate into the solutions of the other components; (2) transferring information between different discretization representations potentially introduces new error; and (3) the adjoint operators associated with the fully coupled system and an operator decomposition version are not generally equal. In addition, the analysis stays within the “single physics paradigm” by only requiring the solution of adjoint problems associated with the individual components. These issues are characteristic of a broad range of operator decomposition discretizations, e.g. [10,13], and generally require extensions to the usual a posteriori analysis techniques. In this paper, we focus attention on analyzing the effects of transferring information between components, which is necessitated by operator decomposition. In order to do so, we consider a “triangular” or one-way coupled system, L1 (x, u1 , Du1 ) = 0 L (x, u1 , Du1 , u2 , Du2 ) = 0 2 L3 (x, u1 , Du1 , u2 , Du2 , u3 , Du3 ) = 0 , x ∈ Ω. (1.5) .. . Ln (x, u1 , Du1 , u2 , Du2 , u3 , Du3 , . . . , un , Dun ) = 0 This system can be solved by a finite sequence of component solutions by considering the n problems for L1 , L2 , · · · , Ln sequentially. Such systems are important in practice, e.g., the thermal actuator (1.1)–(1.3) has this form. In Part II [3], we consider additional sources of error arising from the iterative procedure required when solving a fully-coupled” system via operator decomposition. We capture the essential features of (1.5) in a two component “one-way” coupled system of the form x ∈ Ω, −∇ · a1 ∇u1 + b1 · ∇u1 + c1 u1 = f1 (x), (1.6) −∇ · a2 ∇u2 + b2 · ∇u2 + c2 u2 = f2 (x, u1 , Du1 ), x ∈ Ω, u1 = u2 = 0, x ∈ ∂Ω, where ai , bi , ci , fi are smooth functions on a bounded domain Ω in RN with boundary ∂Ω, and the coupling occurs through f2 . We later generalize to coupling through the coefficients of the elliptic operator for u2 . In Section 2, we illustrate the main idea by applying the analysis to a linear algebraic system. We perform the transfer error analysis in Section 3 and present
4
V. CAREY, D. ESTEP, AND S. TAVENER
computational examples when the corresponding discretizations are “related” in the sense that either both computational meshes are identical, or one mesh is generated by a sequence of mesh refinements on the other mesh. In Section 4, we consider the effect of using distinct discretizations for the two components and analyze the additional errors cause by using projections between the components. Additionally, we discuss the use of Monte Carlo integration to estimate these projection errors. We present the full adaptive algorithm in section 5 which we illustrate with several numerical examples. 2. A linear algebra example. We introduce the notation and ideas in the context of a lower triangular linear system of equations. Let U be an approximate solution of the linear system Au = b. We wish to compute a quantity of interest given by a linear functional (ψ, u). The error, e = u − U is not computable, but we can compute the residual R = b − AU = Ae. Using the solution φ of the corresponding adjoint equation A> φ = ψ, the error representation for a linear functional of the solution is ¡ ¢ ¡ ¢ ¡ ¢ ¡ ¢ ¡ ¢ ¡ ¢ ψ, u − ψ, U = ψ, e = A> φ, e = φ, Ae = φ, R . Now consider the triangular system µ ¶µ ¶ µ ¶ A11 0 u1 b Au = = 1 =b, A21 A22 u2 b2
(2.1)
with approximate solution µ ¶ µ ¶ U1 u1 U= ≈ = u. U2 u2 We estimate the error in a quantity of interest in u2 only, given by the linear functional µ ¶ ¢ ¡ (1) ¢ ¡ (1) 0 . ψ , u = ψ2 , u2 where ψ = (1) ψ2 We employ the superscript (1) since we later pose additional auxiliary adjoint problems. Clearly, estimates on linear functionals of u1 are independent of u2 . The lower triangular structure of A yields A11 u1 = b1 , A22 u2 = b2 − A21 u1 , and the corresponding residuals are R1 = b1 − A11 U1 , R2 = (b2 − A21 U1 ) − A22 U2 . The residual R2 depends upon the solution of the first component, and any attempt to decrease this residual requires a consideration of the accuracy of U1 . The adjoint problem to (2.1) is µ > ¶ Ã (1) ! µ ¶ 0 A11 A> φ1 21 = (1) , (1) 0 A> ψ2 φ2 22
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
and the resulting error representation is ¡ (1) ¡ (1) ¢ ¡ ¢ (1) ψ , e) = ψ2 , e2 = A> 22 φ2 , e2 ¢ ¢ ¡ (1) ¡ (1) = φ2 , A22 u2 − φ2 , A22 U2 ¡ (1) ¢ ¡ (1) ¢ = φ2 , b2 − A21 u1 − φ2 , A22 U2 ¡ (1) ¢ ¡ (1) ¢ = φ2 , b2 − A21 U1 − A22 U2 + φ2 , A21 e1 ¡ (1) ¢ ¡ (1) ¢ = φ2 , R2 + φ2 , A21 e1 .
5
(2.2)
(1)
The first term of the error representation requires only U2 and φ2 . Since the adjoint system is upper triangular and ³ ´−1 (1) (1) φ 2 = A> ψ2 22 ¡ (1) ¢ is independent of the first component, the calculation of φ2 , R2 remains within ¡ (1) ¢ the “single physics paradigm”. The second term φ2 , A21 e1 represents the effect of errors in U1 on the solution U2 . At first glance this term is uncomputable, but we note that it is a linear functional of e1 since ¢ ¢ ¡ ¡ (1) (1) φ2 , A21 e1 = A> 21 φ2 , e1 . We therefore form the adjoint problem for the transfer error, ¶ Ã (2) ! µ (2) ¶ µ > (1) ¶ µ > φ1 A11 A> ψ1 A21 φ2 21 = = . (2) 0 A> 0 0 φ 22 2 (2)
The upper triangular block structure of A> immediately yields φ2 = 0. As noted (2) (2) earlier, error estimates of u1 should be independent of u2 . Thus, A> = 11 φ1 = ψ1 > (1) (2) A21 φ2 , so that once again we can solve for φ in the “single physics paradigm”. Given φ(2) we obtain the secondary error representation ¢ ¡ (2) ¢ ¡ (2) ¢ ¡ > (1) ¢ ¡ > (2) ¢ ¡ (2) (2.3) ψ , e = ψ1 , e1 = A21 φ2 , e1 = A11 φ1 , e1 = φ1 , R1 . Combining the first term of (2.2) with (2.3) yields the complete error representation ¡ (1) ¢ ¡ (1) ¢ ¡ (2) ¢ ψ , e = φ2 , R2 + φ1 , R1 (2.4) which is a sum of the inner products of “single physics” residuals and adjoint solutions computed using the “single physics” paradigm. 3. Analysis of the discretization error. The corresponding weak form of ˜ 1 (Ω) satisfying (1.6) reads: find ui ∈ W 2 ( A1 (u1 , v1 ) = (f1 , v1 ), ˜ 21 (Ω), ∀vi ∈ W (3.1) A2 (u2 , v2 ) = (f2 (x, u1 , Du1 ), v2 ), where A1 (u1 , v1 ) = A1 (u1 , v1 ) ≡ (a1 ∇u1 , ∇v1 ) + (b1 (x) · ∇u1 , v1 ) + (c1 u1 , v1 ), A2 (u2 , v2 ) = A2 (u2 , v2 ) ≡ (a2 ∇u2 , ∇v2 ) + (b2 (x) · ∇u2 , v2 ) + (c2 u2 , v2 ),
6
V. CAREY, D. ESTEP, AND S. TAVENER
˜ m (Ω) is the subspace of W m (Ω) are assumed to be coercive bilinear forms on Ω and W p p with zero trace on ∂Ω. We suppress the “cross” dependence on the other solutions except in a few remarks below. After introducing (conforming) discretizations Sh,i (Ω), we solve the discretized system ( A1 (U1 , χ1 ) = (f1 , χ1 ), ∀χi ∈ Sh,i (Ω). (3.2) A2 (U2 , χ2 ) = (f2 (x, U1 , DU1 ), χ2 ), In general, however, Sh,1 ( Sh,2 (or vice-versa) on Ω, and we may be forced to work with either Π1→2 f2 (U1 ) or more generally with f2 (x, Π1→2 U1 , Π1→2 DU1 ), where ˜ 1 (Ωi ) Πi→j is some projection from Sh,i to Sh,j . If the projection is to Sh,i from W 2 then we simply write the projection as Πi . The resulting discrete system becomes ( A1 (U1 , χ1 ) = (f1 , χ1 ), ∀χi ∈ Sh,i (Ω). (3.3) A2 (U2 , χ2 ) = (f2 (x, Π1→2 U1 , Π1→2 DU1 ), χ2 ), Primary Adjoint Problem We seek the error in a quantity of interest representable by a linear functional of the error e2 , where ui − Ui = ei denotes the pointwise errors. Note that a quantity of interest involving only u1 can be computed without solving for u2 , hence there is no loss of generality. The global adjoint problem, defined relative to the quantity of interest, is ( (1) (1) (1) (1) −∇ · a1 ∇φ1 − div(b1 φ1 ) + c1 φ1 + Lf2 (u1 )φ2 = 0, (1) (1) (1) (1) −∇ · a2 ∇φ2 − div(b2 φ2 ) + c2 φ2 = ψ2 , where Z
1
Lf2 (u1 )(u1 − U1 ) = 0 (1)
∂f2 (u1 s + U1 (1 − s)) ds ∂u1
(1)
is a linearization of f2 and φ1 and φ2 satisfy homogeneous Dirichlet boundary conditions. The corresponding weak formulation is ( (1) (1) A∗1 (φ1 , v1 ) + (Lf2 (u1 )φ2 , v1 ) = 0, ˜ 21 (Ω), ∀vi ∈ W (3.4) (1) (1) A∗2 (φ2 , v2 ) = (ψ2 , v2 ), where ( (1) (1) (1) (1) A∗1 (φ1 , v1 ) = (a1 ∇φ1 , ∇v1 ) − (div(b1 φ1 ), v1 ) + (c1 φ1 , v1 ) (1) (1) (1) (1) A∗2 (φ2 , v2 ) = (a2 ∇φ2 , ∇v2 ) − (div(b2 φ2 ), v2 ) + (c2 φ2 , v2 ).
(3.5)
Using the standard argument, we have the following error representation formula, ¡ (1) ¢ ¡ (1) ¢ ¡ (1) ¢ ¡ ¡ (1) ¢ (1) ¢ ψ , e = ψ2 , e2 = A∗2 φ2 , e2 = f2 (x, u1 , Du1 ), φ2 − A2 U2 , φ2 . (3.6) (1)
Observe that φ1 does not appear in the error representation formula. We define the primary adjoint problem as
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
(1)
7
(1)
˜ 21 (Ω). A∗2 (φ2 , v2 ) = (ψ2 , v2 ), ∀v2 ∈ W Remark 3.1. At first glance, it appears that we only need to solve the second adjoint equation and thus do not need to construct the linearization Lf2 . However, as seen in the linear algebra example, the analysis takes into account the transfer of error from the solution of the first component. Estimating this transferred error uses a nonlinear functional of the error to form the right hand sides in the “transfer adjoint problems”(2.3) and (3.11). We approximate this nonlinear functional using the linearization Lf2 . We evaluate the linearization at the computed solution U , which can be justified by using Taylor’s theorem and assuming that the error u − U is sufficiently small. (1) Adding and subtracting the projection of φ2 onto the primal approximation (1) space (Π2 φ2 ) in (3.6) yields (1)
(1)
(1)
(ψ2 , e2 ) = (f2 (x, u1 , Du1 ), (I − Π2 )φ2 ) − A2 (U2 , (I − Π2 )φ2 ) (1)
(1)
+ (f2 (x, u1 , Du1 ), Π2 φ2 ) − A2 (U2 , Π2 φ2 ). (3.7) To simplify later constructions, we introduce the notion of the weak residual of a solution component, namely Ri (Ui , χ; ν) = (fi (ν), χ) − Ai (Ui , χ; ν), and using this notation write (3.6) as (1)
(ψ (1) , e) = R2 (U2 , φ2 ; u1 ), indicating that this estimate depends on the solution u1 . 3.1. Transfer Error Analysis. The error representation (3.7) is not com¡ (1) ¢ putable since u1 is unknown. We add and subtract f2 (x, U1 , DU1 ), (I − Π2 )φ2 from the error representation formula (3.7) and use the definition of the approximate weak statement (3.2) to obtain ¡ (1) (1) ¢ (1) (ψ2 , e2 ) = f2 (x, U1 , DU1 ), (I − Π2 )φ2 − A2 (U2 , (I − Π2 )φ2 ) ¡ (1) ¢ + f2 (x, u1 , Du1 ) − f2 (x, U1 , DU1 ), φ2 ¢ ¡ ¡ (1) ¢ (1) = R2 U2 , (I − Π2 )φ2 ; U1 + f2 (x, u1 , Du1 ) − f2 (x, U1 , DU1 ), φ2 (3.8) The first term on the right of (3.8) is a traditional dual-weighted residual expression for the error arising from discretization of the second component while the remaining difference represents the transfer error that arises from using an approximation of u1 in defining the coefficients in the equation for u2 . The goal now is to estimate this transfer error and its effect on the quantity of interest. As with the linear algebra example in Section 2, we recognize the transfer error expression as a functional of error in u1 and define (1)
(f2 (x, u1 , Du1 ) − f2 (x, U1 , DU1 ), φ2 ) as a new quantity of interest. Then, we construct a secondary adjoint problem to compute the transfer error. In order to obtain a linear functional when f2 is nonlinear
8
V. CAREY, D. ESTEP, AND S. TAVENER
in u1 , we linearize, f2 (u1 ) ≈ f2 (U1 ) + Df2 (U1 ) × (u1 − U1 ), where Df is the Fr´echet derivative of f2 at U1 . The transfer error term becomes ¡
(1) ¢
Df2 (U1 ) × e1 , φ2
,
(3.9)
which is a linear functional of the error e1 that describes the effect of errors in U1 on the quantity of interest. Note that the Riesz representation theorem guarantees (2) (2) (2) the existence of a ψ1 such that (ψ1 , e1 ) equals (3.9), though ψ1 is not needed to evaluate the functional or compute the corresponding adjoint solution. Transfer Error Adjoint Problem To estimate the new quantity of interest, we define (¡ ¢ ¡ ¢ ¡ (2) ¢ (2) (2) (2) (2) a1 ∇φ1 , ∇v1 − div(b1 φ1 ), v1 + c1 φ1 , v1 ) + (Lf2 (u1 )φ2 , v1 = ψ1 , ¡ ¢ ¡ ¢ ¡ ¢ (2) (2) (2) a2 ∇φ2 , ∇v2 − div(b2 φ2 ), v2 + c2 φ2 , v2 = 0, ˜ 21 (Ω). ∀vi ∈ W The second equation has the trivial solution and the secondary adjoint problem reduces to the “transfer error adjoint problem” ¢ ¡ (2) ¢ ¡ (2) ¢ ¢ ¡ ¡ (2) (2) ˜ 21 (Ω). (3.10) ∀v1 ∈ W a1 ∇φ1 , ∇v1 − (div b1 φ1 , v1 + c1 φ1 , v1 = ψ1 , v1
The transfer error representation formula is given by ¡
¢ (2) (2) (2) ψ1 , e1 = A1∗ (φ1 , e1 ) = A1 (e1 , φ1 ) ¡ (2) (2) ¢ = f1 , (I − Π1 )φ1 − A1 (U1 , (I − Π1 )φ1 ),
(3.11)
where we have used Galerkin orthogonality to introduce the projection of φ onto the discretization space (as f1 does not depend on u). Inserting (3.11) into (3.8) yields ¡ ¢ ¡ (1) ¢ (1) ψ, e = f2 (x, U1 , DU1 ), (I − Π2 )φ2 − A2 (U2 , (I − Π2 )φ2 ) ¢ ¡ (2) (2) + f1 , (I − Π1 )φ1 − A1 (U1 , (I − Π1 )φ1 ). (3.12) or ¡ ¢ (1) (2) ψ, e = R2 (U2 , (I − Π2 )φ2 ; U1 ) + R1 (U1 , (I − Π1 )φ1 ). Remark 3.2. If the model problem includes coupling in the coefficients of the second differential operator, i.e., x ∈ Ω, −∇ · a1 (x)∇u1 + b1 (x) · ∇u1 + c1 (x)u1 = f1 (x), −∇ · a2 (x, u1 )∇u2 + b2 (x, u1 ) · ∇u2 + c2 (x, u1 )u2 = f2 (x, u1 , Du1 ), x ∈ Ω u1 = u2 = 0, x ∈ ∂Ω, (3.13) then the error representation formula for a quantity of interest that depends on u2 alone is ¡ ¢ (1) ψ, e = R2 (U2 , (I − Π2 )φ2 ; u1 ).
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
9
Since this is not computable, we replace each term in the weak residual with the same term evaluated at U1 ; yielding ¡
¢ ¡ (1) (1) ¢ ψ, e = R2 (U2 , (I − Π2 )φ2 ; U1 ) + f2 (u1 ) − f2 (U1 ), φ2 ¡ ¢ ¡ ¢ − (a2 (u1 ) − a2 (U1 ))U2 , φ(1) − (b2 (u1 ) − b2 (U1 )) · ∇U2 , φ(1) ¡ ¢ − (c2 (u1 ) − c2 (U1 ))U2 , φ(1)
We linearize f2 , a2 , b2 , and c2 around U1 to obtain an approximate transfer error term ¡ ¡ ¡ ¡ (1) ¢ (1) ¢ (1) ¢ (1) ¢ Df2 (e1 ), φ2 + Da2 (e1 )∇U2 , ∇φ2 + Db2 (e1 ) · ∇U2 , φ2 + Dc2 (e1 )U2 , φ2 . This is a linear functional on L2 (Ω), which we use as data to define the “transfer” error adjoint problem and derive a corresponding a posteriori error representation. For details on how to compute a quantity of interest that depends on u1 and u2 (so that the choice of linearizations for the coefficients in the equation for u2 enter directly into the “primary” error contribution) see [9]. Remark 3.3. For a “lower triangular” one-way coupled system of N elliptic equations and a quantity of interest based on the N th component, we solve N total “single physics” adjoint problems and construct the error representation ¡
N ¢ X ¡ ¢ ψN , eN = RN −i+1 UN −i+1 , φ(i) ; U . i=1
We then solve a sequence of adjoint problems, as the corresponding linear functional for the ith adjoint problem (i > 1) can be defined recursively (assuming the coupling only occurs through the right-hand side) as ¯ ¶ i−1 µ X DfN +1−j ¯¯ (j) (e ), φ . i Dui ¯ j=1
U
This extends to coupling in all the coefficients as above. 3.2. Numerical examples. The following three numerical examples highlight the features of the analysis and the importance of accounting for the transfer error. In the following computations, we approximately solve all adjoint problems using continuous, piecewise quadratic elements in order to be able to evaluate the interpolants arising from Galerkin orthogonality. We denote these approximate adjoints solutions by Φ and use them in place of φ in the error representation (3.12). For adaptive mesh refinement, we write the estimate as a sum of element contributions and derive a bound by introducing norms. We base the adaptive mesh refinement on the standard optimization approach using the Principle of Equidistribution [6] applied to the bound. We refine elements whose element contribution to the error bound is greater than half a standard deviation from the mean error contribution or refine a fixed fraction of the elements with the greatest element contributions, whichever criterion yields the greater refinement. We do not do any mesh coarsening, smoothing, or edge flips. 3.2.1. Example. This example demonstrates the fact that the transfer error can be significant even if the individual components u1 and u2 are well resolved. We
10
V. CAREY, D. ESTEP, AND S. TAVENER
consider a simple system −∆u1 = sin(4πx) sin(πy), −∆u2 = b · ∇u1 , u = 0, where 2 b= π
µ
25 sin(4πx) sin(πx)
¶
µ ,
f (u) =
(x, y) ∈ Ω, (x, y) ∈ Ω, (x, y) ∈ ∂Ω,
sin(4πx) sin(πy) b · ∇u1
(3.14)
¶ Ω = ([0, 1], [0, 1]).
The quantity of interest is the solution value of u2 at (.25, .25), which we estimate using a smooth delta function approximation with localized support. The corresponding global adjoint problem is (1) (1) −∆φ1 + Lf (u1 )φ2 = 0, (x, y) ∈ Ω (1) (3.15) −∆φ2 = δxreg (x, y) ∈ Ω, ˜ , φ = 0, (x, y) ∈ ∂Ω, where δxreg is a regularized delta function and x ˜ = (.25, .25). Our primary adjoint ˜ problem is (1)
−∆φ2 = δxreg ˜ ,
(x, y) ∈ Ω,
φ = 0, (x, y) ∈ ∂Ω,
The secondary adjoint problem is ( (1) (2) ∆φ1 = ∇ · (bφ2 ), (x, y) ∈ Ω, φ(2) = 0, (x, y) ∈ ∂Ω.
(3.16)
The primal system was solved using identical standard continuous piecewise linear finite element discretizations for u1 and u2 . We plot the results in Fig. 3.1 and show (2) the error contributions in Table 3.1. While the adjoint solution Φ1 in Fig. 3.1(d) is concentrated near the location of the quantity of interest, it has non-trivial spatial structure and the transfer error represents 14% of the total error. Primary Error 0.0042
Transfer Error 0.0006
Table 3.1: Error contributions for Example 3.2.1
3.2.2. Example. This example illustrates the importance of computing the transfer error, since for this problem simply forcing the “primary” error contribution to be small (by refining the second mesh only) does not provide any accuracy in the desired quantity of interest. We reconsider (3.14) but with quantity of interest equal to the average value of u2 over the whole domain. The exact solution has zero average value on Ω. We solve both components of the primary problem on an identical coarse initial mesh, but only adapt and refine the mesh for u2 using the traditional weighted residual, the first “primary” error term in (3.12) while neglecting the second “transfer error” term in (3.12). We show the results in Fig. 3.2 and Table 3.2.
11
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
1 0.5
0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 1
0 −0.5 −1 −1.5 1 0.5 0 0
0.4
0.2
1
0.8
0.6
0.5 0.2
0 0
(a) U1
0.4
0.6
1
0.8
(b) U2
2 0.2 0.15 0.1 0.05 0 −0.05 −0.1 1
1.5 1 0.5 0 1 0.5 0.2
0 0
0.4
0.6
1
0.8
0.5 0.2
0 0
(1)
0.4
0.6
1
0.8
(2)
(c) Φ2
(d) Φ1
Fig. 3.1: Example 3.2.1: Primary and (non-zero) adjoint solutions computed on uniformly fine meshes. The adjoint solutions are largest near the region of the quantity of interest u2 (.25, .25)
0.4 0.2 0.5 0
0
−0.2
−0.5
−0.4
−1
−0.6 1
−1.5 −2 1 0.5 0.5 0
0
0.2
0.4
(a) U1
0.6
0.8
1
0 0
0.2
0.4
0.6
0.8
1
(b) U2
Fig. 3.2: Example 3.2.2: Adaptivity based on the standard discretization error estimate for the “primary” error, ignoring the “transfer” error. Only the mesh for U2 is refined.
12
V. CAREY, D. ESTEP, AND S. TAVENER
Primary Error 0.00005
Transfer Error 0.110
Table 3.2: Error contributions for Example 3.2.2
Ignoring the transfer error and the implied need to refine the first component produces a completely unsuitable adaptive procedure. It is clear from Fig. 3.2 that the average value of the second component is far from zero, and the actual computational value is −0.2245. The estimated transfer error of 0.1 is in fact an underestimate since (2) Φ1 is based on the highly inaccurate solution U1 which is computed on a very coarse mesh. The transfer error dominates the computation and this error cannot be reduced without refining the mesh for u1 . 3.2.3. Example. The third example shows that an “optimal” adaptive mesh for the quantity of interest that depends only on u2 may actually involve a richer discretization of u1 than u2 . We consider the system (3.14) with the quantity of interest equal to the average value of u2 over the whole domain and initial coarse meshes as in the previous example, but use the transfer error contribution to adapt the mesh for u1 and the primary contribution to adapt the mesh for u2 so that the total error is less than 10−4 . The resulting meshes are shown in Fig. 3.3, and illustrate that despite the fact that the quantity of interest only involves u2 , the error inherited from u1 is the most important contribution to consider. In this problem, the strong influence of the transfer error is a result of the dependence of u2 on the gradient of u1 , which a priori has lower order accuracy. Similar behavior could also arise when u2 just depends on u1 .
0.8 0.6 1
0.4
0.5
0.2 0
0 −0.2 −0.5
−0.4 −0.6
−1 1 0.5 0 0
(a) U1
0.2
0.4
0.6
0.8
1
−0.8 1 0.5 0 0
0.2
0.4
0.6
0.8
1
(b) U2
Fig. 3.3: Example 3.2.3: Adaptivity based on the full estimate that accounts for “primary” and “transfer” errors. The quantity of interest U2 (.25, .25) is more sensitive to errors in U1 than U2 .
4. Interpolation error analysis. We use a multiscale discretization for the “fully” adaptive example 3.2.3, i.e., the components u1 and u2 were computed on different meshes, see Figure 3.3. This raises the issue of understanding the effect of translating one component onto the mesh of the other component when performing the
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
13
integration necessary to form the discrete equations. Integration involving functions defined on different meshes these quantities may be complicated, as illustrated in Figure 4.1. Mesh for U1
Mesh for U2
Fig. 4.1: The problem of translation between meshes. Finite element functions on one mesh are generally not smooth on another mesh.
In particular, traditional quadrature formulae based on sets of specific points may not preserve the accuracy required for effective computation because a function defined on a different mesh is generally not sufficiently smooth. For example, the integrand (f2 , χ) is piecewise discontinuous on every element τi of mesh 2 in Example 3.2.3 as b·∇U1 is continuous only within elements of the mesh for U1 . In general, if the meshes are not congruent, the integrand is C 0 at best. Using a “traditional” higher order quadrature rule will not necessarily lead to the expected increase in accuracy as the integrand (f2 , χ) does not have sufficient regularity. Possible solutions include either the determination of local intersections of simplices and/or hexahedra or the construction of a global union mesh. However, both solutions are computationally expensive and the global solution often requires several times more memory than the storage of the two individual meshes, especially for three dimensional problems. 4.1. Projections from mesh 1 to mesh 2. Instead of constructing a union mesh, we use a projection Π1→2 from S1,h to S2,h and solve the discrete system given by (3.3). This introduces additional sources of error. Starting from the error representation formula (3.6), we add and subtract (1)
f2 (x, Π1→2 U1 , Π1→2 DU1 , (I − Π2 )φ2 ), yielding (1)
(ψ (1) , e) = (ψ2 , e2 ) ¡ ¡ (1) ¢ (1) ¢ = f2 (x, Π1→2 U1 , Π1→2 DU1 ), (I − Π2 )φ2 − A2 U2 , (I − Π2 )φ2 ¡ (1) ¢ + f2 (x, u1 , Du1 ) − f2 (x, Π1→2 U1 , Π1→2 DU1 ), φ2 . ¡ (1) ¢ Adding and subtracting f2 (x, U1 , DU1 ), φ2 produces ¡ ¡ (1) ¢ (1) ¢ (ψ (1) , e) = f2 (x, Π1→2 U1 , Π1→2 DU1 ), (I − Π2 )φ2 − A2 U2 , (I − Π2 )φ2 ¡ (1) ¢ + f2 (x, u1 , Du1 ) − f2 (x, U1 , DU1 ), φ2 ¡ (1) ¢ + f2 (x, U1 , DU1 ) − f2 (x, Π1→2 U1 , Π1→2 DU1 ), φ2 .
14
V. CAREY, D. ESTEP, AND S. TAVENER
The first two terms on the right represent the primary discretization error for a functional of the the second component, the third term on the right represents the transfer error (3.11), and the fourth term is a new expression that represents the error from the projection Π1→2 . The projection error can be decomposed as ¡ (1) ¢ Π1→2 f2 (x, U1 , DU1 ) − f2 (x, Π1→2 U1 , Π1→2 DU1 ), φ2 ¡ (1) ¢ + (I − Π1→2 ) f2 (x, U1 , DU1 ), φ2 (4.1) The first inner product in (4.1) can be computed (with some effort) on Ω2,h . However, computing the second term raises the same numerical issues that caused the adoption of the projection Π1→2 in the first place! We handle this term using the Monte Carlo techniques described in Section 4.3. 4.2. Projections from mesh 2 to mesh 1. Complications from the use of projections also arise in computations with the solution of the secondary adjoint problem. (1) The secondary adjoint problem domain is Ω1,h but φ2 is computed naturally on Ω2,h . The new error representation formula for the transfer error becomes ¡ ¡ (1) ¢ (1) ¢ Df2 (U1 ) × e1 , Π2→1 φ2 + Df2 (U1 ) × e1 , (I − Π2→1 )φ2 which is the error contribution arising from the transfer as well as an additional term that is large when the approximation spaces are significantly different. For example, this term is important when the original system is multiscale. The implicit ψ (2) for the transfer error adjoint is now ¡ ¡ ¡ (2) ¢ (1) ¢ (1) ¢ f2 (u1 ) − f2 (U1 ), φ2 = Df2 (U1 ) × e1 , Π2→1 φ2 = ψ1 , e1 ¡ (1) ¢ The additional term Df2 (U1 ) × e1 , (I − Π2→1 )φ2 is a linear functional, so we may define an additional “tertiary” adjoint problem to estimate this quantity. Projection{“Tertiary”} Error Adjoint Problem This problem has the same form as the transfer error adjoint (3.10), but with (3) data ψ1 satisfying ¡ (3) ¢ ¡ (1) ¢ ψ1 , e1 = Df2 (U1 ) × e1 , (I − Π2→1 )φ2 . The resulting error representation formula is ¡ (3) ¢ ¡ ¡ (3) (3) ¢ (3) ¢ ψ1 , e1 = f1 , (I − Π1 )φ1 ) − A1 (U1 , (I − Π1 )φ1 = R1 , (I − Π2→1 )φ1 (4.2) The error representation is therefore (1)
(1)
(2)
(3)
(ψ2 , e2 ) = R2 (U2 , (I − Π2 )φ2 ; U1 ) + R1 (U1 , (I − Π1 )(φ1 + φ1 )) (4.3) ¡ ¡ (1) ¢ (1) ¢ + Π1→2 f2 (U1 ) − f2 (Π1→2 U1 ), φ2 + (I − Π1→2 )f2 (U1 ), φ2 . Remark 4.1. Traditional simplex-based numerical integration methods that interrogates U1 at cubature points can be thought of as projecting the integrand f (U1 )χ2 into a specific polynomial space Pτ defined on each simplex τ of the mesh for U2 and then integrating exactly. We may express this “cubature error” as a projection error and construct a corresponding error representation formula in a similar manner. Cubature error resulting from the fact that integration was not performed on a “union” mesh of two piecewise polynomial spaces may always be viewed as projection error.
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
15
Remark 4.2. In this discussion, we assume that the adjoint problems are solved using approximations spaces that are compatible with the corresponding primal approximation space, e.g., using higher-order Lagrange elements on the same mesh. In practice, different meshes may be used for the primal and adjoint solves. However, this introduces new projection operators between the corresponding approximation spaces as well as the additional terms due to the loss of Galerkin orthogonality. We confine ourselves to merely alluding to the notational complexities and length of the resulting error representation. 4.3. Monte Carlo Integration. Interpolation-based projections suffer from mesh-aliasing difficulties. An extreme example is given in Fig. 4.2. 1
U1 Π1→2U1
0
−1
0
0.2
0.4
0.6
0.8
1
Fig. 4.2: Interpolation errors for two meshes on Ω = [0, 1].
For a more practical example, we construct two quasi-uniform, unstructured meshes, 1 and 2, both of size h on Ω = [0, 1] × [0, 1], and take the piecewise linear interpolantR fI1 of the function R f = sin(20hx) sin(20hy) on mesh 1. We first compute Ie = Ω f dx and I1 = Ω fI1 dx exactly and then construct three different approximations Igauss , IΠ and ISamp as follows. 1. Igauss : Using a 3rd order, 4-point quadrature rule [16] on the triangles of mesh 2 by interpolating fI1 at the corresponding quadrature points. 2. IΠ : Projecting fI1 on to mesh 2 by interpolating fI1 at the nodes of mesh 2 and then using exact integration. 3. ISamp : Performing the integration via a uniform weight quadrature rule using the quadrature points corresponding to the 4-point quadrature rule employed by Igauss . We show the accuracy in Table 4.1. Note that the work for all 3 methods is roughly the same. The smallest of the projection errors |I1 −Igauss |, is larger than the interpolation |Ie − I1 | 0.000187
|I1 − Igauss | 0.000246
|I1 − IΠ | 0.0060
|I1 − ISamp | 0.00041
Table 4.1: Errors in various approximations of Ie . error |Ie − I1 |. The error in |I1 − IΠ | is a factor of 10 larger than |I1 − Igauss | and |I1 − ISamp |, which, for this problem, amounts to a factor of h−1 . Note that the 4-point gauss quadrature rule is only slightly more accurate than the sampling rule ISamp . Motivated by the example, we employ pseudo-random Monte Carlo integration using p random uniformly-distributed sample points on the reference element. The
16
V. CAREY, D. ESTEP, AND S. TAVENER
main difficulty (and computational expense) when integrating on Ω2,h is the evaluation of U1 at each random sample point since this involves locating the point in the appropriate element in Ω1,h . Nominally, this process requires (O(N )) operations per sample point, where N is the number of degrees of freedom for U1 , hence O(M N ) operations for the integration, where M is the number of degrees of freedom for U2 . However, this approach may be greatly accelerated by using a geometric implementation of the assembly and point search algorithms.
τi ∈ Ω2,h pi1 pi2 pi+1 1 u1,h
Fig. 4.3: Monte Carlo integration point search
We illustrate the search algorithm in Fig. 4.3. We generate a random integration point p11 in τ1 ∈ Ω2,h , and determine the containing element of Ω1,h . This could potentially involve a full search of Ω1,h but as this is the initial element, a good starting guess for element location could be provided as an input. Once a matching simplex is found in Ω1,h , the computation is performed and the next integration point p12 is generated. Moreover, the last matching simplex is stored, so the geometric search using edge/face neighbors and barycentric coordinates to guide neighbor selection for the next point is very fast. When the integration is finished, we select the next element to be an edge/face neighbor. Now when we generate p21 , we have a good starting point, namely the last match in p1S which should be “close” to the real element containing p21 . The assembly routine keeps selecting edge neighbors until it has looped over all elements recursively. This algorithm works even with a primitive data structure as long as recursion is employed. If the number of mesh elements is large, however, this may not be practical due to recursion limits. A non-recursive algorithm could lead to termination before all element contributions for the mesh were calculated, as the next element returned by the search could have all edge/face neighbors whose element contributions had already been calculated. The algorithm would have to “restart” from an element that has not been computed. On quasi-uniform meshes with no fine scale features in the geometry, the number of “restarts” also grows logarithmically with the number of elements. Of course, with a more sophisticated data structure, either octree based or, for example, a mesh where the elements had been ordered by the use of a space filling curve, the need for restarting would be eliminated. When the meshes for Ω1,h and Ω2,h are both quasi-uniform on Ω, the number of elements tested in Ω is bounded by some h-independent constant for each integration point. Obviously this is not the case for general adapted or anisotropic meshes, but in practice the number of searches grows at most logarithmically with the numbers of degrees of freedom in u1 . The convergence of this Monte Carlo integration scheme follows from standard results (see [11]) as the integrand can always be defined as the
17
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
sum of integrals of continuous functions on individual simplices of the union mesh of Ω1,h and Ω2,h . 4.4. Numerical examples. We demonstrate the significance of the projection errors with two examples. 4.4.1. Example. The first example illustrates how the projection error can influence a typical computation. We consider a system defined by (3.14), with two randomly-generated initial meshes for u1 and u2 . The initial mesh for u1 is finer than for u2 in order to reduce the transfer error. The quantity of interest in this computation is the average value of u2 . We show the results in Fig. 4.4 and Table 4.2. We use a local projector Π1→2,τ given by interpolation at the Gauss points (third order 3 point simplex rule) of simplices τ in Sh,2 . Use of this projector would integrate (U1 , U2 ) exactly if the meshes were identical. The solution using this projector is given by Fig. 4.4(b). This is compared against a 16 point Monte Carlo computation illustrated by Fig. 4.4(c).
2 1.5
1
1 0.5
0.5 0
0
−0.5 −0.5
1 0.8
−1 1
0.6 1
0.4
0.8
0.5
0.6
0.2
0.4 0
0.2
0
0
(a) U1
0.2
0
0.4
0.6
0.8
1
(b) U2 computed with Π1→2,τ
0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 1
0.5
0
0
0.2
0.4
0.6
0.8
1
(c) U2 computed with high-sample Monte Carlo integration
Fig. 4.4: Example 4.4.1: The role of projection errors in non-aligned meshes. Note that the magnitude and oscillation of U2 computed with Π1→2,τ , shown in (b) are incorrect (a fine scale U2 is given by Fig. 3.1).
18
V. CAREY, D. ESTEP, AND S. TAVENER
Primary Error 0.003533
Transfer Error 0.021589
Projection Error 0.007908
Table 4.2: Example 4.4.1: Error contributions for computation shown in Fig. 4.4.
2 2 As discussed in Section 4.2, projection from Sh,2 to Sh,1 can also lead to significant inaccuracies in computing the transfer error, necessitating the computation of a tertiary adjoint problem (4.2). 2 2 to Sh,1 can 4.4.2. Example. As discussed in Section 4.2, projection from Sh,2 also lead to significant inaccuracies in computing the transfer error, necessitating the computation of a tertiary adjoint problem (4.2). This example shows that computations with significant differences in mesh scale can contribute significantly to the error. We again use the system in Example 3.2.1 with quantity of interest the point value at (.15, 15), starting with a coarse identical initial mesh for u1 and u2 , but refining only the mesh for u2 . There is no projection error as Sh,2 ⊆ Sh,1 . However, when we compute the transfer error, we ignore the fact that a natural choice of decomposition for the computation is integration over the simplices of Sh,2 . Instead, we use (1) interpolation of φ2 at the quadrature points at the simplices of Sh,1 . To compute (I − Π), we employ the actual nesting of the two meshes to perform an accurate (up (3) to quadrature error on the fine scale mesh) computation of φ1 . We show the results in Table 4.3 and Fig. 4.5.
Primary Error 0.000713
Transfer Error 0.0905
Projection Error 0
Tertiary Error 0.0325
Table 4.3: Example 4.4.2: Error contributions for the computation shown in Fig. 4.5.
1 0.5 0 −0.5 1 −1 0
0.5 0.2
0.4
0.6
0.8
1 0
(3)
Fig. 4.5: Example 4.4.2: Tertiary adjoint solution Φ1 which estimates the projection error in computing the transfer error
5. An adaptive algorithm for the operator decomposition - finite element method. An adaptive algorithm that takes into account all the possible
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
19
sources of error is given below. while (the total error is less than TOL) do Compute U1 using standard integration. Compute U2 using 16 point M.C. integration for the coupling term. (1) Compute Φ2 using standard integration (2) Compute Φ1 for given adjoint data using 16 point M.C. integration if (the sum of two error contributions is greater than TOL) then Refine both meshes based on the primary error contributions for U2 and the transfer error contributions for U1 . else Compute the projection error by comparing with a 64 point M.C. integration. (3) Compute Φ1 . if (the total error is greater than TOL) then Refine both meshes based on the primary and projection error contributions for U2 and the transfer and tertiary error contributions for U1 . end if end if end while The algorithm drives the primary and transfer error contributions to within a specified error tolerance, and then checks for projection error by using 64 point Monte Carlo integration as an approximation to the identity operator I in (4.3) and attempts to correct the projection error by refinement as well. Any projector could be substi(2) tuted for the M.C. integration used in computing U2 and Φ1 . We select the use of 16 sample points for the Monte Carlo integration based on our experience from a series of numerical experiments where different functions were interpolated on a quasi uniform mesh, and integrated. This interpolant was then integrated using Monte Carlo with 2N sample points per simplex on a different quasiuniform mesh (with the same approximate h); N = 4 gave the best tradeoff between speed and accuracy. 5.1. Examples. We describe two applications of the algorithm to one way coupled systems using different meshes for each solution component. In both examples, we start with identical coarse initial meshes (quasi-uniform with h ≈ .125) and adapt each mesh until both the primary and transfer error formulas are less than 10−4 . We control projection error using Monte Carlo integration. 5.1.1. Example. In the first example, we approximate the value of u2 at (.25, .25), where (u1 , u2 ) solves 2 −∆u1 = 64π sin 4π(x − .75 + |x − .75|) sin 4π(y − .75 + |y − .75|), −∆u2 = u1 , u1 = u2 = 0, with Ω = ([0, 1], [0, 1]). The corresponding adjoint problem is ( (1) −∆φ2 = δreg (x0 ), (x, y) ∈ Ω, (1) φ2 = 0, (x, y) ∈ ∂Ω,
(x, y) ∈ Ω, (x, y) ∈ Ω, (x, y) ∈ ∂Ω,
20
V. CAREY, D. ESTEP, AND S. TAVENER
with x0 = (.25, .25). The transfer error adjoint problem is ( (2) (1) −∆φ1 = φ2 , (2) φ1 = 0,
(x, y) ∈ Ω, (x, y) ∈ ∂Ω.
1
1
0.5 0.5 0 1 0.8 0
0.6 0.4 0.2
−0.5 1
0.8
0.6
0.4
0.2
0 1
0.5
0
0
0
0.2
0.4
0.6
0.8
1
(a) U1 computed using a uniformly fine mesh, (b) U2 for the quantity of interest u2 (0.25, 0.25) axes rotated
0.8 0.6 0.4 0.2
1 0.8
0
0.6
−0.2
0.4
−0.4
0.2
−0.6
0
0
−0.4 1
1
0.2
−0.2 0.4
0.5
0.6 0.8
0.6
0.4
0.8 0.2
0
1
0 0
0.2
0.4
0.6
0.8
1
(c) U1 for the quantity of interest u2 (0.25, 0.25), (d) U1 for the quantity of interest u2 (0.9, 0.9) axes rotated
Fig. 5.1: Example 5.1.1: Example of computational efficiency: U1 may be computed on a coarse discretization yet U2 may be determined with sufficient accuracy.
The accurate computation of the quantity of interest u2 (0.25, 0.25) does not require the fine scale features of u1 near (0.75, 0.75) to be resolved. However, a quantity of interest equal to the value of u2 at (0.9, 0.9), near the localized features of u1 requires better resolution of the details of u1 . The adapted solutions U1 for both quantities of interest and given in Fig. 5.1(c) and Fig. 5.1(d) respectively. 5.1.2. Example. We now consider an example where convection in component u1 creates the need for refinement in u1 remote from the the goal-oriented refinement
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
in u2 .
3 −100kx−x0 k2 , x ∈ Ω, −∆u1 − b · ∇u1 = 10 e 3 −100kx−x1 k2 −∆u2 = 10 e u1 , x ∈ Ω, u1 = u2 = 0, x ∈ ∂Ω
21
(5.1)
where b = (100 40)> , x0 = (.75, .75), x1 = (.1, .5), and the quantity of interest is the point value u2 (x2 ), x2 = (.2, .5). The corresponding adjoint problem for the primary error contribution is ( (1) −∆φ2 = δreg (x2 ), x ∈ Ω, (1) φ2 = 0, x ∈ ∂Ω, while the corresponding transfer error adjoint problem is ( 2 (1) (2) (2) −∆φ1 + b · ∇φ1 = 103 e−100kx−x1 k φ2 , x ∈ Ω, (2) φ1 = 0, x ∈ ∂Ω. (2)
The adjoint solution Φ1 in Fig. 5.2(c) shows the influence of the convection term in the equation for u1 . When the quantity of interest is a value of u2 in the convective region of influence of the localized source term in the equation for u1 , the solution for u1 is resolved “upstream” of the location of the quantity of interest as shown in Fig. 5.2(a). When the quantity of interest is a value of u2 away from the convective region of influence of the localized source term in the equation for u1 , the adjoint solution (2) φ1 has an similar structure to that shown in Fig. 5.2(c) but has much smaller magnitude. The resulting mesh for U1 need not even be detailed enough to eliminate the numerical oscillation (from not satisfying the corresponding P´eclet mesh condition). This situation is illustrated by Fig. 5.2(d), where the choice of quantity of interest is u2 at (0.15, 0.15). 6. Conclusion. In this paper, we perform an a posteriori error analysis of a multiscale operator decomposition finite element method for the solution of a system of one-way coupled elliptic problems. The analysis specifically accounts for the effects arising from multiscale operator decomposition, including the issues: (1) errors in the solution of each component propagate into the solutions of the other components; and (2) transferring information between different representations potentially introduces new error. We estimate the various sources of errors by defining auxiliary adjoint problems whose data are related to errors in the information passed between components. Through a series of examples, we demonstrate the importance of accounting for the contributions to the error arising from multiscale operator decomposition. We also devise an adaptive discretization strategy based on the error estimates that specifically controls the effects arising from operator decomposition. Finally, we demonstrate the usefulness of Monte Carlo integration methods for dealing with mismatch between discretizations of different components. We extend this analysis to a “fully-coupled” system in the form of (1.4) in Part II of this paper [3]. We address the important issue that the adjoint operator associated with the fully coupled system and an operator decomposition solution are not generally equal. This difference requires additional strategies for error control. We consider the use of non-interpolatory projectors based on averaging to reduce both transfer and projection error in [4].
22
V. CAREY, D. ESTEP, AND S. TAVENER
1.5 5 1 4 0.5 3 0 2 −0.5 1
1 0 1 0.5 0 0
0.2
0.4
0.8
0.6
1
0.5 0 0
(a) U1 for quantity of interest u2 (.2, .5)
0.2
0.4
0.6
0.8
1
(b) U2 for quantity of interest u2 (.2, .5)
0.5 1.5
0
1
−0.5 1
0.5
0.8
0
0.6 −0.5 1
0.4
1
0.2
0.8
0.5
0.6 0.4
0
0 (2)
(c) Φ1
0.2
0.4
0.6
0.8
1
for quantity of interest u2 (.2, .5)
0
0.2 0
(d) U1 for quantity of interest u2 (.15, .15)
Fig. 5.2: The role of convection in Example 5.1.2. Note that altering the location of the quantity of interest alters the density and location of the resulting adapted meshes (same adaptive criteria).
REFERENCES [1] Wolfgang Bangerth and Rolf Rannacher, Adaptive Finite Element Methods for Differential Equations, Birkhauser Verlag, 2003. [2] R. Becker and R. Rannacher, An optimal control approach to a posteriori error estimation in finite element methods, Acta Numerica, (2001), pp. 1–102. [3] Varis Carey, Donald Estep, and Simon Tavener, A Posteriori analysis and adaptive error control for operator decomposition methods for elliptic systems II: Fully coupled systems, In preparation, (2007). [4] , Averaging based projections in operator decomposition methods for elliptic systems, In preparation, (2007). [5] K. Eriksson, D. Estep, P. Hansbo, and C. Johnson, Introduction to adaptive methods for differential equations, in Acta numerica, 1995, Acta Numer., Cambridge Univ. Press, Cambridge, 1995, pp. 105–158. [6] , Computational Differential Equations, Cambridge University Press, Cambridge, 1996. [7] D. Estep, A posteriori error bounds and global error control for approximation of ordinary differential equations, SIAM J. Numer. Anal., 32 (1995), pp. 1–48. [8] D. Estep, M. Holst, and M. Larson, Generalized Green’s functions and the effective domain of influence, SIAM J. Sci. Comput., 26 (2005), pp. 1314–1339. [9] D. Estep, M. G. Larson, and R. D. Williams, Estimating the error of numerical solutions of
MULTISCALE OPERATOR DECOMPOSITION FOR ELLIPTIC SYSTEMS
23
systems of reaction-diffusion equations, Mem. Amer. Math. Soc., 146 (2000), pp. viii+109. [10] D. Estep, S. Tavener, and T. Wildey, A posteriori analysis and improved accuracy for an operator decomposition solution of a conjugate heat transfer problem. SINUM, in revision, 2006. [11] George S. Fishman, Monte Carlo, Springer Series in Operations Research, Springer-Verlag, New York, 1996. Concepts, algorithms, and applications. ¨ li, Adjoint methods for PDEs: A posteriori error analysis and postpro[12] M. Giles and E. Su cessing by duality, Acta Numerica, (2002), pp. 145–236. [13] V. Ginting, D. Estep, J. Shadid, and S. Tavener, An a posteriori analysis of operator splitting. to appear in SINUM, 2006. [14] G.I. Marchuk, On the theory of the splitting-up method, in Proceedings of the Second Symposium on Numerical Solution of Partial Differential Equations, SVNSPADE, Academic Press, New York, 1970, pp. 469–500. [15] , Splitting and alternating direction methods, in Handbook of Numerical Analysis, vol. I, North-Holland, New York, 1990, pp. 197–462. [16] Gilbert Strang and George J. Fix, An Analysis of the Finite Element Method, Prentice-Hall Inc., Englewood Cliffs, N. J., 1973. Prentice-Hall Series in Automatic Computation.
