A POSTERIORI ANALYSIS AND IMPROVED ... - Semantic Scholar
A POSTERIORI ANALYSIS AND IMPROVED ACCURACY FOR AN OPERATOR DECOMPOSITION SOLUTION OF A CONJUGATE HEAT TRANSFER PROBLEM D. ESTEP∗ , S. TAVENER† , AND T. WILDEY‡ Abstract. We consider the accuracy of an operator decomposition finite element method for a conjugate heat transfer problem consisting of two materials coupled through a common boundary. We derive accurate a posteriori error estimates that account for the transfer of error between components of the operator decomposition method as well as the differences between the adjoints of the full problem and the discrete iterative system. We use these estimates to guide adaptive mesh refinement. In addition, we address a loss of order of convergence that results from the decomposition, and show that the approximation order of convergence is limited by the accuracy of the transferred gradient information. We employ a boundary flux recovery method to regain the expected order of accuracy in an efficient manner. Key words. a posteriori error analysis, adaptive mesh refinement, adjoint problem, boundary flux method, conjugate heat transfer, domain decomposition, finite element method, generalized Green’s function, goal oriented error estimates, operator decomposition, residual, transfer error AMS subject classifications. 65N15, 65N30, 65N50
1. Introduction. In this paper, we consider the solution of a conjugate heat transfer problem by an operator decomposition approach. The goal is to compute a functional of the temperature of a body composed of two distinct components that share a common boundary or interface. The two components may have different conductivities and be subject to different heat sources and boundary conditions. The model consists of a system of boundary value problems for two elliptic equations coupled through boundary conditions posed on the common interface. Our interest is in situations where there is a small number of interfaces between the components. Operator decomposition provides an alternative strategy to solving the fully coupled system. We compute the temperature in each component individually subject to the conditions on the common interface that are obtained from a numerical solution of the other component. Obtaining a solution of the full system by operator decomposition therefore requires a (nominally) infinite iteration in which updated boundary conditions are passed between components. Operator decomposition is a widely used technique for solving multi-physics, multi-scale problems. The general approach is to decompose the problem into components involving simpler physics over a relatively limited range of scales, and then to seek the solution of the entire system through an iterative procedure involving solutions of the individual components. This approach is appealing because there is ∗ 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), the National Aeronautics and Space Administration (NNG04GH63G), the National Science Foundation (DMS-0107832, DGE0221595003, MSPA-CSE-0434354) and the Sandia Corporation (PO299784) † Department of Mathematics, Colorado State University, Fort Collins, CO 80523 ([email protected]). S. Taveners’s work is supported in part by the Department of Energy (DE-FG02-04ER25620) and the National Science Foundation (xxx) ‡ Department of Mathematics, Colorado State University, Fort Collins, CO 80523 ([email protected]). T. Wildey’s work is supported in part by the Department of Energy (DE-FG02-04ER25620), the National Science Foundation (DMS-0107832) and the Sandia Corporation (PO299784)
1
2
D. ESTEP, S. TAVENER, AND T. WILDEY
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. In the case of conjugate heat transfer, operator decomposition can be viewed as domain decomposition that allows for very different discretizations in each component. However, 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. In the case of the conjugate heat transfer, the operator decomposition causes a loss in the approximation order of convergence. In this paper, we perform an a posteriori error analysis of a finite-element implementation of the operator decomposition technique and obtain accurate error estimates that are used to guide an adaptive discretization strategy. Our approach is based on the standard techniques using variational analysis, residuals and the generalized Green’s function solution to an adjoint problem [3, 6, 14, 15, 12, 20], which we modify to account for several new features arising from the operator decomposition. These include the following. • Numerical errors in the solution of each component are propagated through the boundary condition applied to the other component. • Numerical errors in the solution of each the components at one step of the iterative procedure are propagated to the next step. • The adjoint of the full problem and the adjoint of the iterative procedure are distinct operators whose differences must be recognized when seeking to obtain accurate error estimates. These kinds of effects are characteristic of operator decomposition discretizations, e.g. [8, 19], and generally require extensions of the usual a posteriori analysis techniques. In addition to obtaining accurate estimates, we seek to improve the accuracy of the operator decomposition method in an efficient way. In particular, we adapt the “boundary element flux” technique developed by Wheeler [25] and Carey [16, 7] to compute normal derivatives on a boundary, and show that this can be used to improve accuracy, and in particular, restore the order of convergence that is lost due to the operator decomposition. The operator decomposition method we consider is closely related to the nonoverlapping domain decomposition methods as described in, for example, [23, 28, 24, 29]. The major concern of these authors is however the convergence of different iterative methods and the design of effective pre-conditioners or relaxation parameters, while the accuracy of the final finite element approximation is usually not addressed, despite the fact that some of the numerical results indicate a loss of accuracy. To our knowledge, this paper is the first to address the loss of accuracy. Our results would extend to this particular type of non-overlapping domain decomposition. In §2, we introduce the conjugate heat transfer problem and provide some notation. We describe the iterative operator decomposition finite element method and some modifications in §3, as well as the boundary flux method used to compute gradients on the common interface. We perform two a posteriori error analyses in §4, using first the adjoint to the fully coupled problem and then the adjoint to the iterative scheme, and present numerical results to highlight the differences between the two. We end the section with a brief discussion of adaptive mesh refinement. In §5, we carry out an analysis which identifies the transferred gradient information as being responsible for the loss of order and show that using the boundary flux method to compute the gradient information restores the order of convergence. Our conclusions
3
CONJUGATE HEAT TRANSFER
are presented in §6. 2. The model for conjugate heat transfer. Let Ω1 and Ω2 be polygonal domains in R2 or R3 with boundaries ∂Ω1 and ∂Ω2 intersecting along an interface Γ = ∂Ω1 ∩ ∂Ω2 . We consider a system of second order linear elliptic problems, where the components are coupled through boundary conditions imposed on Γ, L1 u1 = f1 , u (1 = 0, u =u , 1 2 A1 ∂n u1 = A2 ∂n u2 , L2 u2 = f2 , u2 = 0,
x ∈ Ω1 , x ∈ ∂Ω1 \Γ, x ∈ Γ,
(2.1)
x ∈ Ω2 , x ∈ ∂Ω2 \Γ,
where for i = 1, 2, Li ui = −∇ · (Ai ∇ui ) + ci ui , Ai ≥ Ai,0 > 0, ci , fi are sufficiently smooth functions and ∂n denotes the unit outward normal derivative to ∂Ω1 . The results of this paper extend easily to general elliptic operators and general Dirichlet, Neumann, and Robin boundary conditions on the boundaries in the complement of the interface. We let L2 (Ωi ) denote the space of square integrable functions on Ωi with inner product (·, ·)Ωi and norm k · kΩi , but use (·, ·) = (·, ·)Ωi when the domain is clear. We use H s (Ωi ) to denote the Sobolev space with real index s associated with the norm k · kΩi ,s and seminorm | · |Ωi ,s . We also use the subspaces H01 (Ωi ) = 1 v ∈ H (Ωi ), v = 0 on ∂Ωi \Γ . The weak formulation of (2.1) seeks ui ∈ H01 (Ωi ) such that u1 = u2 on Γ and 2 X i=1
ai (ui , vi ) =
2 X
(fi , vi ),
(2.2)
i=1
R for all vi ∈ H01 (Ωi ) with ai (ui , v) = Ωi (Ai ∇ui · ∇v + ci ui v) dx for i = 1, 2. Assuming that each ai (·, ·) is coercive, (2.2) admits a unique weak solution in H01 (Ωi ) [1, 5, 9]. 3. An iterative operator decomposition method. To compute a numerical (0) solution of (2.1), we consider iterative operator decomposition methods. Let u2 be an initial guess for the Dirichlet data along the interface. We solve (k) x ∈ Ω1 , L1 u1 = f1 , (k) (3.1) u1 = 0, x ∈ ∂Ω1 \Γ, (k) (k−1) u1 = u2 , x ∈ Γ, followed by
(k) L2 u2 = f2 , (k) u2 = 0, (k) (k) A1 ∂n u1 = A2 ∂n u2 ,
x ∈ Ω2 , x ∈ ∂Ω2 \Γ, x ∈ Γ.
(3.2)
(k)
We iterate until (hopefully) a convergence criteria is met, e.g., ku2 smaller than a prescribed tolerance.
(k−1)
− u2
kΓ is
4
D. ESTEP, S. TAVENER, AND T. WILDEY
3.1. Finite element discretization. We let Ti,h be a triangulation of Ωi into elements K where the length of the longest edge is hK and h = maxK∈Ti,h hK . We assume that each triangulation is locally quasi-uniform and Ωi = ∪K∈Ti,h K. However, the triangulations on either side of Γ are not assumed to be aligned. We use the piecewise polynomial spaces S1 = v continuous on Ω1 , v ∈ P 1 (K) for all K ∈ T1,h , S2 = v continuous on Ω2 , v ∈ P 1 (K) for all K ∈ T2,h , and the associated spaces
S1,0 = {v ∈ S1 | v = 0, x ∈ ∂Ω1 } , S2,0 = {v ∈ S2 | v = 0, x ∈ ∂Ω2 \Γ} , where P 1 (K) denotes the space of linear polynomials on an element K. We let πi be a projection into Si as well as the projection into Si along the interface Γ. (k) For the finite element approximation, we compute U1 ∈ S1 satisfying (k) a1 (U1 , v) = (f1 , v)Ω1 , for all v ∈ S1,0 , (k) (3.3) x ∈ ∂Ω1 \Γ, U1 = 0, (k) (k−1) U1 = π1 U2 , x ∈ Γ, (k)
followed by U2 ∈ S2 such that ( (k) (k) a2 (U2 , v) = (f2 , v)Ω2 − (A1 ∂n U1 , v)Γ , (k) U2 = 0,
for all v ∈ S2,0 , x ∈ ∂Ω2 \Γ.
(3.4)
3.2. Relaxed iterations. Unfortunately, the simple iterative (3.3)-(3.4) may not converge. In particular, the convergence depends on the values of A1 and A2 along the interface and the geometry of each region [18, 23, 28]. As an alternative, we consider two “relaxed” iteration schemes. (1) Relaxed Dirichlet values We choose α ∈ [0, 1] and update the Dirichlet interface values with (k)
U1
(k−1)
= αU1
(k−1)
+ (1 − α)π1 U2
.
(3.5)
Optimal values of α can be found in [23, 28], but α = 1/2 works well in most situations. (2) Relaxed Neumann values (k−1) (k) We chose β ∈ [0, 1], set Nβ = βA2 ∂n U2 + (1 − β)A1 ∂n U1 , and seek (k) U2 ∈ S2 such that (k)
a2 (U2 , v) = (f2 , v)Ω2 − (Nβ , v)Γ , for all v ∈ S2,0 ,
(3.6)
A proper choice of β reduces the number of iterations. 3.3. Flux correction. It turns out that the operator decomposition results in a loss of order arising from the relatively low order of accuracy in the derivative of the finite element solution, which is O(h) when computed with a standard O(h2 ) method. We show below that this reduces the order of the overall approximation to first order.
CONJUGATE HEAT TRANSFER
5
To mitigate this effect, we use a post-processing technique developed by Wheeler [25] and Carey [16, 7] to compute more accurate boundary flux information. We define the set of elements in T1,h that intersect the boundary by Γ T1,h = {K ∈ T1,h | K ∩ ∂Ω 6= ∅},
and the corresponding space Γ Wh = {v ∈ P 1 (K) with K ∈ T1,h , v(ηi ) = 0 if ηi ∈ / ∂Ω},
where {ηi } denotes the nodes of element K, so the degrees of freedom correspond to the nodes on the boundary. We seek σ (k) ∈ Wh satisfying (k)
−(σ (k) , v)Γ = (f1 , v)Ω1 − a1 (U1 , v), for all v ∈ Wh ,
(3.7)
(k)
where U1 solves (3.3). Green’s identity implies that σ (k) gives an approximation to the normal flux on the boundary which is relatively inexpensive to compute. In general, the accuracy of the flux approximation depends on the regularity of an associated Green’s function [17, 27], but in many cases this post-processed flux is O(h2 ) rather than the standard O(h) for the normal flux of a piecewise linear finite element approximation. It turns out that a fortunate cancellation of errors makes the accuracy of this recovered flux unimportant for our purposes. We stress that σ (k) is not assumed to be continuous. In fact, if the domain has a corner on the boundary, the normal derivative is, in general, discontinuous due to the jump in the normal vector. When Dirichlet conditions are given on each boundary segment we implement the method described in [7] and allow two degrees of freedom to account for the discontinuity. For the application to more general boundary conditions, such as a Neumann or Robin condition on one segment, or an interface condition on both, we refer the reader to [26]. Now, as an alternative to (3.4), we may solve ( (k) a2 (U2 , v) = (f2 , v)Ω2 − (σ (k) , v)Γ , for all v ∈ S2,0 , (3.8) (k) U2 = 0, x ∈ ∂Ω2 \Γ. Other possible approaches to mitigating the loss of order include approximating the boundary flux using a gradient recovery techniques such as the ZienkiewiczZhu (ZZ) patch recovery technique [30] or the polynomial preserving recovery (PPR) method [22], or using higher order polynomials near the interface to improve the accuracy of the finite element flux [2, 20]. We had mixed computational success using these other approaches. 4. A posteriori error analysis. To estimate the error of the operator decomposition finite element approximation, we apply a posteriori techniques based on variational analysis and the adjoint problem. In this case, however, the adjoint for the fully coupled original problem differs significantly from the adjoint problem associated with the discretization of the decomposed system. The motivation to use operator decomposition suggests that the adjoint of the full problem will be unavailable in practice and therefore the decomposed adjoint will be used to compute error estimates. The difference between the adjoint problems can be decreased by increasing the computational work used to solve the decomposed adjoint, however in practice, we must consider the issue that the two adjoints lead to significantly different error representations.
6
D. ESTEP, S. TAVENER, AND T. WILDEY
In order to understand the differences, we derive a posteriori error estimates using both adjoints. Both analyses begin in the same way. We wish to estimate the difference between the exact solution to the full problem and the numerical approximation to the solution of the iterative procedure, i.e. E (k) = u − U (k) . We decompose the error, E (k) = u − U (k) = u − u(k) + u(k) − U (k) = c(k) + e(k) .
The first component c(k) measures the difference between the exact solution to the full problem and the exact solution to the iterative problem. If the iterative method converges, then c(k) → 0. The second component e(k) measures the numerical error in solving the iterative problem. It is helpful to remark that in the conventional approach, the adjoint is defined in such a way as to make the formal bilinear identity hold [21, 13], i.e., to make the boundary terms arising from integration by parts equal to zero. In the analysis below, we find it convenient to consider the values passed between components as quantities of interest defined on the common boundary and we have to alter the definition of the adjoint problem accordingly. 4.1. The adjoint to the fully coupled problem. The adjoint boundary value problem for the quantity of interest (ψ, u) = (ψ1 , u1 )+(ψ2 , u2 ) for the coupled problem (2.1) is ∗ x ∈ Ω1 , L1 φ1 = ψ1 , φ x ∈ ∂Ω1 \Γ, (1 = 0, φ =φ , 1 2 (4.1) x ∈ Γ, A ∂ φ 1 n 1 = A2 ∂n φ2 , L∗2 φ2 = ψ2 , x ∈ Ω2 , φ2 = 0, x ∈ ∂Ω2 \Γ, where L∗i φi = −∇ · (Ai ∇φi ) + ci φi . We use a∗i (·, ·) to represent the weak form of the adjoint operator. We solve (4.1) numerically by using an iterative operator decomposition approach as for the forward problem. The iterations are completely independent of the forward iterations. We can derive an error representation formula for the basic scheme (3.1)-(3.2), a weighted relaxation technique (3.5) or (3.6), or when using the post-processed flux as (k) in (3.8). In the discussion below, we use θh to denote the numerical flux passed at the k th iteration from Ω1 to Ω2 . To begin, we multiply the system (4.1) by (ψ1 , ψ2 )T and apply the divergence theorem, noting that u1 = u2 and A1 ∂n φ1 = A2 ∂n φ2 along Γ, to obtain (ψ1 , e1 ) + (ψ2 , e2 ) = a1 (e1 , φ1 ) + a2 (e2 , φ2 ) (k)
(k)
+ (U1 , A1 ∂n φ1 )Γ − (U2 , A2 ∂n φ2 )Γ . Observe that the test space S1,0 consists of functions that are zero along the interface, while in general, φ1 is not zero along Γ. This means that the projection of φ1 into S1,0 cannot be the interpolant. We define a new projection π10 : H 2 → S1,0 such that for any node xi ( π1 φ(xi ), xi ∈ / Γ, π10 φ(xi ) = (4.2) 0, xi ∈ Γ.
7
CONJUGATE HEAT TRANSFER
We also observe that (k)
a2 (e2 , π2 φ2 ) = −(A1 ∂n u1 , π2 φ2 )Γ + (θh , π2 φ2 )Γ , and a1 (u1 , φ1 − π10 φ1 ) = (f1 , φ1 − π10 φ1 ) + (A1 ∂n u1 , φ1 )Γ , a2 (u2 , φ2 − π2 φ2 ) = (f2 , φ2 − π2 φ2 ) − (A1 ∂n u1 , φ2 − π2 φ2 )Γ , since the adjoint solutions are not zero along Γ. Using the projection (4.2) in the Galerkin orthogonality relation and the equalities above, we have (k)
(ψ1 , e1 ) + (ψ2 , e2 ) = (f1 , φ1 − π10 φ1 ) − a1 (U1 , φ1 − π10 φ1 ) (k)
+ (f2 , φ2 − π2 φ2 ) − a2 (U2 , φ2 − π2 φ2 ) (k)
+ (U1
(k)
(k)
− U2 , A1 ∂n φ1 )Γ + (θh , π2 φ2 )Γ .
Next, we define π∂ φ1 = π1 φ1 − π10 φ1 which is nonzero only near the interface due to the definition of π10 φ1 . Substituting π10 φ1 = π1 φ1 − π∂ φ1 gives (k)
(ψ1 , e1 ) + (ψ2 , e2 ) = (f1 , φ1 − π1 φ1 ) − a1 (U1 , φ1 − π1 φ1 ) (k)
+ (f2 , φ2 − π2 φ2 ) − a2 (U2 , φ2 − π2 φ2 ) (k)
+ (f1 , π∂ φ1 ) − a1 (U1 , π∂ φ1 ) (k)
+ (U1
(k)
(k)
− U2 , A1 ∂n φ1 )Γ + (θh , π2 φ2 )Γ
Finally (3.7) implies that the recovered boundary flux, σ (k) satisfies (k)
−(σ (k) , π∂ φ1 )Γ = (f1 , π∂ φ1 )Ω1 − a1 (U1 , π∂ φ1 ), while π∂ φ1 = π1 φ1 along Γ. We conclude (k) (k) Theorem 4.1. The errors e1 = u1 − U1 and e2 = u2 − U2 satisfy (k)
(ψ1 , e1 ) + (ψ2 , e2 ) = (f1 , φ1 − π1 φ1 ) − a1 (U1 , φ1 − π1 φ1 ) (k)
+ (f2 , φ2 − π2 φ2 ) − a2 (U2 , φ2 − π2 φ2 ) (k)
+ (U1
(k)
− U2 , A1 ∂n φ1 )Γ
(k)
+ (θh , π2 φ2 )Γ − (σ (k) , π1 φ1 )Γ
(4.3)
The error has been decomposed into two typical discretization components, an iterative component, and a component reflecting the error arising from the transfer of derivative information. The choice of iterative method does not affect the discretization component of the error, but it does influence both the iterative and the transfer components along Γ. We illustrate this for a few common iterative schemes: (k) (k−1) • Suppose U1 = π1 U2 , then (k)
(U1
(k)
− U2 , A1 ∂n φ1 )Γ (k−1)
= (π1 U2
(k−1)
− U2
(k)
, A1 ∂n φ1 )Γ + (U2
(k−1)
− U2
which represents a projection error and an iteration error.
, A1 ∂n φ1 )Γ ,
8
D. ESTEP, S. TAVENER, AND T. WILDEY (k)
• Suppose U1 (k)
(U1
(k−1)
= αU1
(k−1)
+ (1 − α)π1 U2
(k)
− U2 , A1 ∂n φ1 )Γ =
k−2 X
, then (k−i−1)
αi−1 (π1 (U1
(k−i)
− U2
), A1 ∂n φ1 )Γ
i=1
(k−1)
+ (U2
(k)
− U2 , A1 ∂n φ1 )Γ
(k−1)
+ (π1 U2
(k−1)
− U2
, A1 ∂n φ1 )Γ ,
which represents a series of iteration errors and a projection error. Notice that since α < 1, the effect of the iteration error from previous iterations diminishes due to the increasing power on α. (k) (k) • Suppose we set θh = A1 ∂n U1 , then (k)
(k)
(θh , π2 φ2 )Γ −(σ (k) , π1 φ1 )Γ = (A1 ∂n U1 −σ (k) , π2 φ2 )Γ +(σ (k) , π2 φ2 −π1 φ1 )Γ , which represents a transfer error and a projection error. (k) • Suppose we set θh = σ (k) , then (k)
(θh , π2 φ2 )Γ − (σ (k) , π1 φ1 )Γ = (σ (k) , π2 φ2 − π1 φ1 )Γ , which represents only a projection error with no transfer error. 4.2. The adjoint to the iterative scheme. Next, we derive an error representation using the natural adjoint for the iterative system (3.1)-(3.2), where we set (k) (k−1) (k) (k) U1 = π1 U2 and A2 ∂n U2 = A1 ∂n U1 . This adjoint avoids the formulation of a globally coupled adjoint, and reads (k) ∗ (k) L1 φ1 = ψ1 , x ∈ Ω1 , (k) φ1 = 0, x ∈ ∂Ω1 \Γ, (k) (k) φ1 = φ2 , x ∈ Γ, (k) ∗ (k) L2 φ2 = ψ2 , (k) φ2 = 0, (k) (k+1) A2 ∂n φ2 = A1 ∂n φ1 ,
x ∈ Ω2 , x ∈ ∂Ω2 \Γ, x ∈ Γ,
(N +1)
(4.4)
(4.5)
for k = N, . . . , 1 with A1 ∂n φ1 = 0. Note that the adjoint system is defined “backwards”. In the reasonable case in which we seek information from the final (1) (1) (N −1) iterate, the data for the adjoint problem would be ψ1 = ψ2 = . . . = ψ1 = (N −1) (N ) (N ) ψ2 = 0, with ψ1 and ψ2 chosen appropriately. We derive the error representation formula for this system by observing N X
k=1
N X (k) (k) (k) (k) (k) (k) (k) (k) (L∗1 φ1 , e1 ) + (L∗2 φ2 , e2 ) , (ψ1 , e1 ) + (ψ2 , e2 ) = k=1
and apply the steps used to derive (4.3). This provides the following theorem.
9
CONJUGATE HEAT TRANSFER (1)
(1)
(N )
(N )
Theorem 4.2. The errors e1 , e2 , . . . , e1 , e2 N X
k=1
satisfy
(k) (k) (k) (k) (ψ1 , e1 ) + (ψ2 , e2 )
=
N X
(k)
(k)
(k)
(k)
(k)
(f1 , φ1 − π1 φ1 ) − a1 (U1 , φ1 − π1 φ1 )
k=1
(k)
(k)
(k)
(k)
(k)
+ (f2 , φ2 − π2 φ2 ) − a2 (U2 , φ2 − π2 φ2 ) +
(k) (U1
−
(k−1) (k) U2 , A1 ∂n φ1 )Γ
+
(k) (k) (A1 ∂n U1 , π1 φ1 )Γ
− (σ
(k)
(k) , π2 φ2 )Γ
(4.6) It is possible to define adjoints for more complicated iterative methods, but it (k) (k−1) may be difficult. For example, relaxation techniques couple U1 to both U2 and (k−1) U1 which affects how the adjoint variables are coupled. In addition, using the post-processed numerical flux given in section 3.2 requires defining the adjoint of the post-processing procedure. 4.3. Numerical results. We illustrate the a posteriori estimates in §4.1 and §4.2. Example 4.1. We triangulate the domains Ω1 = [−0.25, 0.25] × [−0.25, 0.25] and Ω2 = ([−1, 1] × [−1, 1]) \Ω1 independently, see Fig. 4.1, and consider the elliptic
Fig. 4.1. Triangulations of Ω1 and Ω2 that do not match along the interface.
interface problem, −∇ · (∇u1 ) = f1 , ( u1 = u2 , ∂n u 1 = ∂n u 2 , −∇ · (∇u2 ) = f2 , u2 = 0,
x ∈ Ω1 , x ∈ Γ,
(4.7)
x ∈ Ω2 , x ∈ ∂Ω2 \Γ,
where f1 (x, y) and f2 (x, y) are chosen so the true solutions are u2 = sin(4πx) sin(4πy) (k) (k−1) and u1 = 3u2 . We solve this using (3.3)-(3.4) with U1 = π1 U2 and iterate until (k) (k−1) −5 kU1 − U1 kΓ < 10 . We consider several quantities of interest
10
D. ESTEP, S. TAVENER, AND T. WILDEY
1. 2. 3. 4.
the global average value, the average value over Ω1 , the average value over Ω2 , the value of u2 at the point (0, 0) computed using the approximate delta 2 2 function δˆ0 (x, y) = 400 π exp(−400x − 400y ). First we use an independent iterative scheme to solve the adjoint of the fully (k) coupled problem using the relaxation parameter α = 1/2 and iterating until kΦ1 − (k−1) (k) Φ1 kΓ < 10−6 , where Φ1 represents the approximation of the adjoint solution at th the k iteration. We display the error estimates for each of the quantities of interest in Table 4.1. Next we solve the adjoint to the operator decomposition method and give these results in Table 4.2. Note that the number of iterations, as well as the values of the relaxation parameters, are determined by the number of iterations used in the solution of the forward operator decomposition method.
1. 2. 3. 4.
ψ1 4 4 0 δˆ0
ψ2 4/15 0 4/15 0
True Error −0.20 −0.15 −0.046 −0.15
F.C. Adjoint −0.20 −0.15 −0.046 −0.15
Effect. Ratio 0.99 0.99 1.00 0.99
Table 4.1 Error estimates and effectivity ratios using the adjoint for the fully coupled system.
1. 2. 3. 4.
ψ1 4 4 0 δˆ0
ψ2 4/15 0 4/15 0
True Error −0.20 −0.15 −0.046 −0.15
O.D. Adjoint −0.20 −0.15 −0.046 −0.15
Effect. Ratio 0.99 0.99 1.00 0.99
Table 4.2 Error estimates and effectivity ratios using the adjoint for the operator decomposition method.
Example 4.2. Next we consider neighboring domains Ω1 = [0, 1]×[0, 1] and Ω2 = ([1, 2] × [0, 1]) with independent triangulations, and the elliptic interface problem −∇ · (A1 ∇u1 ) = f1 , u (1 = 0, u =u , 1 2 A ∂ u 1 n 1 = A2 ∂n u2 , −∇ · (A2 ∇u2 ) = f2 , u2 = 0,
x ∈ Ω1 , x ∈ ∂Ω1 \Γ, x ∈ Γ,
(4.8)
x ∈ Ω2 , x ∈ ∂Ω2 \Γ,
where A1 = 1, A2 = 3, and f1 (x, y) and f2 (x, y) are chosen so the true solutions are u2 = sin(2πx) sin(2πy) and u1 = 3u2 . The quantity of interest is the average value over Ω1 ∪ Ω2 , so ψ1 = ψ2 = 1. This experiment demonstrates that the error representation formulas (4.3) and (4.6) may give different estimates if relatively few iterations are used in the operator
11
CONJUGATE HEAT TRANSFER 1.2 1.6
Adjoint for Fully Coupled Problem Adjoint for Operator Decomposition Method
Size of adjoints
Effectivity Ratios
1.1
1
0.9
1.2
0.8
0.4 0.8 2
4
6
8 10 Iterations
12
14
0
16
0
5
10
15
Iteration
Fig. 4.2. On the left, a comparison of the effectivity ratios using the two adjoints for a given number of iterations used for the operator decomposition method for the forward problem. On the (k) (k) right, we plot kφ1 k+kφ2 k to show the decay of influence of errors that occur in previous iterations.
decomposition method for the forward problem. We set a fixed number of iterations in the operator decomposition method for (4.8) and use (3.3) and (3.4) to approximate the solution. We solve the fully coupled adjoint iteratively using sufficient iterations for convergence. However, the number of iterations for the adjoint of the iterative system is determined by the number of iterations for the forward problem. In Fig. 4.2, we plot the number of iterations versus the effectivity ratios for the two estimates. We observe that the adjoint for the operator decomposition method does not produce accurate estimates until a sufficient number of iterations for the forward problem have been carried out. This implies that the iterative component of the error, c(k) = u − u(k) , dominates the estimate. If the iterative method converges, we expect that the effect of the initial guess and the first few iterations would diminish. This is reflected in Fig. 4.2 where we (k) (k) plot the L2 norm of the iterative adjoint φ1 and φ2 over Ω1 and Ω2 respectively at each iteration. We observe that the norm decays rapidly as k decreases, which indicates that the norm of the fixed point operator, and hence the norm of the adjoint fixed point operator, is small. Thus, the influence of errors in previous iterations are rapidly damped and have little influence on the current error. This implies that we can compute accurate estimates using the truncated series (N )
(ψ1
(N )
(N )
(N )
, e1 ) + (ψ2 , e2 ) N X (k) (k) (k) (k) (k) (f1 , φ1 − π1 φ1 ) − a1 (U1 , φ1 − π1 φ1 ) = k=M
(k)
(k)
(k)
(k)
(k)
+ (f2 , φ2 − π2 φ2 ) − a2 (U2 , φ2 − π2 φ2 ) (k)
+ (U1
(k−1)
− U2
(k)
(k)
(k)
(k)
, A1 ∂n φ1 )Γ + (A1 ∂n U1 , π1 φ1 )Γ − (σ (k) , π2 φ2 )Γ
(4.9) where 1 ≤ M ≤ N , in place of the full series (4.6). This significantly reduces the computational cost associated with computing a series of adjoint problems. In Fig. 4.3 we plot the effectivity ratio for the global average value where we truncate the error representation and compute only the last N terms.
12
D. ESTEP, S. TAVENER, AND T. WILDEY
Effectivity Ratio
1.15 1.05 0.95 0.85 0.75 0.65 0
2
4
6 8 10 Number of Terms
12
14
Fig. 4.3. A comparison of the effectivity ratios computed using a truncated error representation.
4.4. Adaptive refinement. We use the a posteriori error estimate as the basis for adaptivity by employing the standard “optimization framework” after writing the estimate as a sum of element contributions and introducing norms [10, 11, 4, 3]. An element K is marked for refinement when the local error indicator is larger than a local tolerance, usually the global tolerance divided by the current number of elements. Example 4.3. We demonstrate the adaptive procedure on the model problem (4.7) with the quantity of interest equal to the value of u2 at the point (1.75, 0.25), which we approximate by choosing ψ1 = 0 and ψ2 = 400/π exp(−400(x − 1.75)2 − 400(y − 0.25)2 ). Fig. 4.4 shows the results after 3 refinement steps. We solve (3.3)(3.4) and observe that refinement occurs within Ω1 along Γ which reflects the impact of numerical errors in the normal derivative on the quantity of interest.
Fig. 4.4. Adaptive mesh for the quantity of interest equal to the value of u2 at the point (1.75, 0.25).
5. An analysis of the loss of order. In practice, the iterative technique (3.3) and (3.4) is occasionally observed to result in O(h) convergence rather than the O(h2 ) convergence that is obtained when solving the fully coupled problem. This loss of order results from passing the normal derivative of the finite element approximation, which is only O(h). Passing the recovered boundary flux restores the full order of convergence. We use the adjoint for the fully coupled problem to derive a posteriori (k) (k−1) error bounds for the iterative approximations in the case where U1 = π1 U2 . The case where the Dirichlet values are updated using a relaxation technique can be analyzed using the same approach and gives identical results. Numerical examples are provided at the end of the section.
CONJUGATE HEAT TRANSFER
13
5.1. L2 error bounds. Let u1 ∈ H 1+α1 (Ω1 ) and u2 ∈ H 1+α2 (Ω2 ) be the solu(k) (k) tions to (2.1) with 0 ≤ α1 , α2 ≤ 1, and U1 and U2 be the solutions of (3.3) and th 1+α1 (Ω1 ) and φ2 ∈ H 1+α2 (Ω2 ) and (3.4) respectively at the k iteration. Let φ1 ∈ H pose the adjoint problem (4.1) with ψ1 = e1 /ke1 kΩ1 and ψ2 = e2 /ke2 kΩ2 . Starting with (4.3), integration by parts over each element K gives ke1 kΩ1 + ke2 kΩ2 = I1 + I2 + I3 + I4 , where I1 =
X
K∈T1,h
1 (k) (k) (f1 − L1 U1 , φ1 − π1 φ1 )K + ([A1 ∂n U1 ], φ1 − π1 φ1 )∂K 2 X
+
K∈T2,h
1 (k) (k) (f2 − L2 U2 , φ2 − π2 φ2 )K + ([A2 ∂n U2 ], φ2 − π2 φ2 )∂K 2
(k)
(k)
I2 =(A2 ∂n U2 , φ2 − π2 φ2 )Γ − (A1 ∂n U1 , φ1 − π1 φ1 )Γ , (k)
(k)
− U2 )Γ ,
I3 =(A1 ∂n φ1 , U1 (k)
I4 =(A1 ∂n U1 , π2 φ2 )Γ − (σ (k) , π1 φ1 )Γ , where [Ai ∂n Ui ] denoting the jump in the normal derivative across an element edge. The first term I1 is a standard a posteriori error bound for elliptic problems and is not affected by non-matching triangulations along the interface or by transfer error. The second term I2 is similar to the jump terms along element edges in I1 and is the expected jump term when the triangulations align along the interface. The third term I3 represents the jump in the Dirichlet values across the interface. Finally, the fourth term I4 represents the difference between the flux passed from Ω1 to Ω2 and the flux obtained via the boundary-flux correction technique. We first construct Lemmas 5.1 to 5.4 below to bound I1 to I4 individually. In each of these Lemmas we first provide the general bound when the triangulations do not match across the boundary and then show the simplification that arises for matching triangulations. We then combine these four lemmas into Theorems 5.5 and 5.7 which give error bounds for the basic iteration (3.3) and (3.4), and when using flux correction (3.8) respectively. These two theorems describe the general result when the triangulations do not match across the boundary while the simplification given matching triangulations is provided as a Corollary. Lemma 5.1. (Bound on I1 ) I1 ≤
X
K∈T1,h
+
X
! 1 Ch1+α |φ1 |K,1+α1 K · 1 Ch1+α |φ1 |K,1+α1 K ! (k) 2 kf2 − L2 U2 kK Ch1+α |φ2 |K,1+α2 K . · −1/2 (k) 2 1 Ch1+α |φ2 |K,1+α2 [A2 ∂n U2 ]k∂K K 2 khK (k)
kf1 − L1 U1 kK −1/2 (k) 1 [A1 ∂n U1 ]k∂K 2 khK
K∈T2,h
Proof. Apply standard arguments using the Cauchy-Schwarz inequality and a trace inequality.
14
D. ESTEP, S. TAVENER, AND T. WILDEY
If u1 ∈ H 2 (Ω1 ), we can also provide a bound for the jump terms. Ordinarily, we would add and subtract (A1 ∂n u1 , φ1 − π1 φ1 )∂K to estimate the jump term. However, (k) this causes an iteration error (A1 ∂n u1 − A1 ∂n u1 , φ1 − π1 φ1 )Γ , which may be quite (k) complicated. Therefore, we let u ˜1 solve
(k)
(k) u1 , v) = (f1 , v)Ω1 , a1 (˜ (k) u ˜1 = 0, (k) (k−1) u ˜1 = π1 U2 ,
for all v ∈ H01 (Ω1 ), x ∈ ∂Ω\Γ, x ∈ Γ.
(5.1)
To be precise, u ˜1 is the exact weak solution over Ω1 given the data from the previous iteration. Standard finite element theory can be applied to estimate (k)
(k)
k˜ u1 − U1 kΩ1 ,1 ≤ Ch1 kf1 kΩ1 . (k)
We add and subtract (A1 ∂n u ˜1 , φ1 − π1 φ1 )∂K and use the Cauchy-Schwarz and trace inequalities to obtain, (k)
1/2
k[A1 ∂n U1 ]k∂K ≤ ChK kf1 kΩ1 . A similar argument holds if u2 ∈ H 2 (Ω2 ). Lemma 5.2. (Bound on I2 ) If the triangulations T1,h and T2,h do not match along the interface, then (k) −1/2 2 |φ2 |1+α2 + kA1 ∂n U1 kΓ · (kπ1 φ1 − π2 φ1 kΓ ) I2 ≤ kh2 [A∂n U (k) ]kΓ · Ch1+α 2
where
kπ1 φ1 − π2 φ1 kΓ ≤ kπ1 φ1 − φ1 kΓ + kφ1 − π2 φ1 kΓ 1/2+α1
≤ C1 h1
1/2+α2
|φ1 |Ω1 ,1+α1 + C2 h2
|φ2 |Ω2 ,1+α2 .
If the triangulations T1,h and T2,h match along the interface, then −1/2 2 |φ2 |1+α2 . I2 ≤ kh2 [A∂n U (k) ]kΓ · Ch1+α 2 (k)
where [A∂n U (k) ] = A2 ∂n U2
(k)
− A1 ∂n U1
and h2 = maxK∈T2,h hK . (k)
Proof. We add and subtract (A1 ∂n U1 , φ2 − π2 φ2 )Γ and use φ1 = φ2 on Γ to write (k)
I2 = (A2 ∂n U2
(k)
(k)
− A1 ∂n U1 , φ2 − π2 φ2 )Γ + (A1 ∂n U1 , π1 φ1 − π2 φ1 )Γ .
Observe π1 φ1 = π2 φ2 if the triangulations match along the interface. We apply the Cauchy-Schwarz and trace inequalities to complete the proof. (k)
If u2 ∈ H 2 (Ω2 ), we may define u ˜2
(k)
similarly to u ˜1 1/2
k[A∂n U ]kΓ ≤ Ch2 kf2 kΩ1 .
in (5.1) and obtain
15
CONJUGATE HEAT TRANSFER
Lemma 5.3. (Bound on I3 ) If the triangulations T1,h and T2,h do not match along the interface, then (k−1) (k) I3 ≤ (kA1 ∂n φ1 kΓ ) · kU2 − U2 kΓ (k−1) (k−1) + (kA1 ∂n φ1 kΓ ) · kπ1 U2 − U2 kΓ where
(k−1)
kπ1 U2
(k−1)
− U2
1/2+α1
kΓ ≤ C1 h1
1/2+α2
|ˆ u|Ω1 ,1+α1 + C2 h2
|ˆ u|Ω2 ,1+α2 .
If the triangulations T1,h and T2,h match along the interface, then (k−1) (k) I3 ≤ (kA1 ∂n φ1 kΓ ) · kU2 − U2 kΓ . (k)
Proof. First, observe that U1 (k−1) tracting (A1 ∂n φ1 , U2 )Γ to get (k−1)
I3 = (A1 ∂n φ1 , π1 U2 (k−1)
(k−1)
= π1 U2 (k−1)
− U2
and rewrite I3 by adding and sub(k−1)
)Γ + (A1 ∂n φ1 , U2
(k)
− U2 )Γ .
(k−1)
To bound kπ1 U2 − U2 kΓ we introduce a new function u ˆ such that u ˆ ∈ (k−1) (k−1) H 1+α1 (Ω1 ) and u ˆ ∈ H 1+α2 (Ω2 ), as well as π1 u ˆ = π1 U2 and π2 u ˆ = U2 . This yields (k−1)
kπ1 U2
(k−1)
− U2
(k−1)
1/2+α1
kΓ ≤ C1 h1
1/2+α2
|ˆ u|Ω1 ,1+α1 + C2 h2
|ˆ u|Ω2 ,1+α2 .
(k−1)
Observe that π1 U2 = U2 if the triangulations match along the interface. We apply the Cauchy-Schwarz inequality to complete the proof. Lemma 5.4. (Bound on I4 ) If the triangulations T1,h and T2,h do not match along the interface, then (k) (k) (k) I4 ≤ kA1 ∂n U1 − A1 ∂n u ˜1 kΓ + kA1 ∂n u ˜1 − σ (k) kΓ · (kπ2 φ2 kΓ ) + kσ (k) kΓ · (kπ2 φ2 − π1 φ1 kΓ )
where
1/2+α1
kπ2 φ2 − π1 φ1 kΓ ≤ C1 h1
1/2+α2
|φ1 |Ω1 ,1+α1 + C2 h2
|φ2 |Ω2 ,1+α2 .
If the triangulations T1,h and T2,h match along the interface, then (k) (k) (k) I4 ≤ kA1 ∂n U1 − A1 ∂n u ˜1 kΓ + kA1 ∂n u ˜1 − σ (k) kΓ · (kπ2 φ2 kΓ ) . (k)
where u ˜1
is defined by (5.1).
Proof. We add and subtract (σ (k) , π2 φ2 )Γ and use φ1 = φ2 along Γ to get (k)
I4 = (A1 ∂n U1
− σ (k) , π2 φ2 )Γ + (σ (k) , π2 φ2 − π1 φ2 )Γ .
16
D. ESTEP, S. TAVENER, AND T. WILDEY
Observe π1 φ1 = π2 φ2 if the triangulations match along the interface. We add and (k) subtract (A1 ∂n u ˜1 , π2 φ2 )Γ to the first term and apply the Cauchy-Schwarz inequality to complete the proof.
In practice, the error in the normal derivative is typically the same accuracy as 1 the H 1 error, namely O(hα 1 ). However, an application of the trace theorem only α1 /2 proves O(h1 ) accuracy. This is not an important issue, however, since we intend to use the fact that this term is less accurate than the others, and therefore pollutes the L2 error. We assume the error in the normal derivative can be bounded, (k)
kA1 ∂n U1
(k)
−u ˜1 kΓ ≤ Chβ1 kf1 kΩ1 ,
for α1 /2 ≤ β ≤ α1 . The error in the recovered boundary flux can be bounded (k)
1 kf1 kΩ1 , kA1 ∂n u ˜1 − σ (k) kΓ ≤ CS1 h1+α 1
where S1 is a stability factor defined by an associated Green’s function [17, 27].
Theorem 5.5. Assume the triangulations T1,h and T2,h do not match along (k) (k) the interface Γ. If U1 and U2 solve (3.3) and (3.4) respectively, then the errors (k) (k) e1 = u1 − U1 and e2 = u2 − U2 satisfy
ke1 kΩ1 + ke2 kΩ2 ≤
X
K∈T1,h
+
X
K∈T2,h
! 1 Ch1+α |φ1 |K,1+α1 K · 1 Ch1+α |φ1 |K,1+α1 K ! (k) 2 kf2 − L2 U2 kK Ch1+α |φ2 |K,1+α2 K · −1/2 (k) 2 1 Ch1+α |φ2 |K,1+α2 [A2 ∂n U2 ]k∂K K 2 khK −1/2 2 |φ2 |1+α2 + kh2 [A∂n U (k) ]kΓ · Ch1+α 2 (k−1) (k) + (kA1 ∂n φ1 kΓ ) · kU2 − U2 kΓ 1 + Chβ1 kf1 kΩ1 + CS1 h1+α kf1 kΩ1 · (kπ2 φ2 kΓ ) 1 (k)
kf1 − L1 U1 kK −1/2 (k) 1 [A1 ∂n U1 ]k∂K 2 khK
+C1 h1
1/2+α1
u|Ω1 ,1+α1 ) ((S2 + S4 )|φ1 |Ω1 ,1+α1 + S3 |ˆ
1/2+α2 +C2 h2
u|Ω2 ,1+α2 ) ((S2 + S4 )|φ2 |Ω2 ,1+α2 + S3 |ˆ (k)
with α1 /2 ≤ β ≤ α1 , S1 is a stability factor independent of h, S2 = kA1 ∂n U1 kΓ , S3 = kA1 ∂n φ1 kΓ , and S4 = kσ (k) kΓ .
Corollary 5.6. Assume the triangulations T1,h and T2,h match along the (k) (k) interface Γ. If U1 and U2 solve (3.3) and (3.4) respectively, then the errors
CONJUGATE HEAT TRANSFER (k)
e1 = u1 − U1
(k)
and e2 = u2 − U2
ke1 kΩ1 + ke2 kΩ2 ≤
satisfy
! 1 Ch1+α |φ1 |K,1+α1 K · 1 Ch1+α |φ1 |K,1+α1 K ! (k) 2 kf2 − L2 U2 kK Ch1+α |φ2 |K,1+α2 K · −1/2 (k) 2 1 Ch1+α |φ2 |K,1+α2 [A2 ∂n U2 ]k∂K K 2 khK −1/2 2 |φ2 |1+α2 + kh2 [A∂n U (k) ]kΓ · Ch1+α 2 (k−1) (k) + (kA1 ∂n φ1 kΓ ) · kU2 − U2 kΓ 1 kf1 kΩ1 · (kπ2 φ2 kΓ ) , + Chβ1 kf1 kΩ1 + CS1 h1+α 1 (k)
X
K∈T1,h
+
17
X
K∈T2,h
kf1 − L1 U1 kK −1/2 (k) 1 kh [A1 ∂n U1 ]k∂K K 2
with α1 /2 ≤ β ≤ α1 and S1 is a stability factor independent of h1 . It is clear that the term containing hβ1 decreases at a slower rate than the other terms. In fact, if we assume u1 ∈ H 2 (Ω1 ) and u2 ∈ H 2 (Ω2 ) then the other terms are O(h21 ) or O(h22 ), while the finite element flux is O(h1 ) at best. Suppose that instead we solve (3.8) rather than (3.4), i.e., we pass σ (k) instead of the finite element flux. This changes the fourth term in the error representation formula to I4 = (σ (k) , π2 φ2 )Γ − (σ (k) , π1 φ1 )Γ .
Theorem 5.7. Assume the triangulations T1,h and T2,h do not match along (k) (k) the interface Γ. If U1 and U2 solve (3.3) and (3.8) respectively, then the errors (k) (k) e1 = u1 − U1 and e2 = u2 − U2 satisfy ke1 kΩ1 + ke2 kΩ2 ≤
K∈T1,h
+
X
! 1 Ch1+α |φ1 |K,1+α1 K · 1 Ch1+α |φ1 |K,1+α1 K ! (k) 2 Ch1+α |φ2 |K,1+α2 kf2 − L2 U2 kK K · −1/2 (k) 2 1 Ch1+α |φ2 |K,1+α2 [A2 ∂n U2 ]k∂K K 2 khK −1/2 2 |φ2 |1+α2 + kh2 [A∂n U (k) ]kΓ · Ch1+α 2 (k−1) (k) + (kA1 ∂n φ1 kΓ ) · kU2 − U2 kΓ (k)
X
K∈T2,h
kf1 − L1 U1 kK −1/2 (k) 1 kh [A1 ∂n U1 ]k∂K K 2
1/2+α1
+C1 h1
1/2+α2 +C1 h2
u|Ω1 ,1+α1 ) ((S2 + S4 )|φ1 |Ω1 ,1+α1 + S3 |ˆ
u|Ω2 ,1+α2 ) , ((S2 + S4 )|φ2 |Ω2 ,1+α2 + S3 |ˆ
(k)
where S2 = kA1 ∂n U1 kΓ , S3 = kA1 ∂n φ1 kΓ , and S4 = kσ (k) kΓ . Corollary 5.8. Assume the triangulations T1,h and T2,h match along the (k) (k) interface Γ. If U1 and U2 solve (3.3) and (3.8) respectively, then the errors
18
D. ESTEP, S. TAVENER, AND T. WILDEY (k)
e1 = u1 − U1
(k)
and e2 = u2 − U2
ke1 kΩ1 + ke2 kΩ2 ≤
X
K∈T1,h
+
X
K∈T2,h
satisfy
! 1 Ch1+α |φ1 |K,1+α1 K · 1 Ch1+α |φ1 |K,1+α1 K ! (k) 2 Ch1+α |φ2 |K,1+α2 kf2 − L2 U2 kK K · −1/2 (k) 2 1 Ch1+α |φ2 |K,1+α2 [A2 ∂n U2 ]k∂K K 2 khK −1/2 2 |φ2 |1+α2 + kh2 [A∂n U (k) ]kΓ · Ch1+α 2 (k−1) (k) + (kA1 ∂n φ1 kΓ ) · kU2 − U2 kΓ . (k)
kf1 − L1 U1 kK −1/2 (k) 1 kh [A1 ∂n U1 ]k∂K K 2
Comparing Theorem 5.7 with Theorem 5.5 and Corollary 5.8 with Corollary 5.6, we see that the terms involving hβ have dropped out and the optimal order of convergence has been restored. 5.2. Numerical results. Example 5.1. Let Ω1 = [0, 1] × [0, 1] and the triangulations match along Γ = {(x, y)| x = interface problem −∇ ( · (∇u1 ) = f1 , u =u , 1 2 ∂ u = ∂n u 2 , n 1 −∇ · (∇u2 ) = f2 ,
Ω2 = [1, 2] × [0, 1] and assume that 1, 0 ≤ y ≤ 1}. Consider the elliptic x ∈ Ω1 , x ∈ Γ,
(5.2)
x ∈ Ω2 ,
where the data f1 (x, y) and f2 (x, y) are chosen so the true solutions are u1 = u2 = sin(πx/2) sin(2πy). The purpose of this experiment is to show how the asymptotic O(hβ1 ) convergence rate may be masked by apparently minor changes to the problem, and therefore may be difficult to observe in all situations. We change the boundary conditions along the top boundary, solving (5.2) with the boundary conditions u1 = 0, along x = 0 and y = 0, ∂ u = 2π sin(πx/2), along y = 1, n 1 u along x = 2 and y = 0, 2 = 0, ∂n u2 = 2π sin(πx/2), along y = 1. We set α = 1/2 and iterate until the convergence criteria is met. For comparison, we also solve (5.2) using a fully coupled method on the same mesh in COMSOL. We see, in Figure 5.1, that the fully coupled solution converges quadratically, while the iterative solution converges linearly. Next, we solve (5.2) with Dirichlet boundary conditions ( u1 = 0, x ∈ ∂Ω1 \Γ, u2 = 0, x ∈ ∂Ω2 \Γ, and again we also solve the fully coupled problem for comparison. In Figure 5.2, we observe that the operator decomposition solution initially converges quadratically, but eventually the O(hβ1 ) term dominates and the convergence is ultimately linear.
19
CONJUGATE HEAT TRANSFER
Fully Coupled Solution Operator Decomposition
Fully Coupled Solution Operator Decomposition
−3
−4
log( L2error)
log( L2 error)
−3
−5 −6
−4 −5 −6
−7 1.5
2 2.5 3 log(sqrt(degrees of freedom))
−7 1.5
3.5
2 2.5 3 log(sqrt(degrees of freedom))
3.5
Fig. 5.1. Comparison of L2 error in the fully coupled approximation and the operator decomposition approximation over Ω1 (left) and Ω2 (right) with mixed boundary conditions on ∂Ωi \Γ for i = 1, 2. −2
−2 Fully Coupled Solution Operator Decomposition
Fully Coupled Solution Operator Decomposition
−4 log(L2error)
log(L2 error)
−4 −6 −8 −10
−6 −8 −10
2
3
4
5
6
log(sqrt(degrees of freedom))
2
3 4 5 log(sqrt(degrees of freedom))
6
Fig. 5.2. Comparison of L2 error in the fully coupled approximation and the iterative approximation over Ω1 (left) and Ω2 (right) with Dirichlet boundary conditions on ∂Ωi \Γ for i = 1, 2.
Example 5.2. In the final example we assume the triangulations do not match along the interface and consider (5.2) with the mixed boundary conditions u1 = 0, A1 ∂n u1 A ∂ u 1 n 1 u = 0, 2 A2 ∂n u2 A2 ∂n u2
= 2π sin(πx/2), = −2π sin(πx/2), = 2π sin(πx/2), = −2π sin(πx/2),
along along along along along along
x = 0, y = 1, y = 0, x = 2, y = 1, y = 0. (k)
First, we solve the problem iteratively by passing the finite element flux, A1 ∂n U1 . Next, we use the boundary flux method to compute and pass σ (k) . In Figure 5.3, we compare the L2 errors over Ω1 ∪ Ω2 on a series of meshes. We see that the approximations using the finite element flux converge linearly, while the approximations using the boundary flux converge quadratically. 6. Conclusion. We have conducted an a posteriori error analysis of a finite element implementation of a number of operator decomposition strategies for a canonical
20
D. ESTEP, S. TAVENER, AND T. WILDEY
Using Finite Element Flux Using Boundary Flux
−5
2
log(L error)
−4
−6
−7 2.8
3
3.2 3.4 3.6 3.8 4 log(sqrt(degrees of freedom))
4.2
4.4
Fig. 5.3. Comparison of L2 error in the iterative approximation using the finite element flux and the boundary flux recovery method over Ω1 ∪ Ω2 with mixed boundary conditions on ∂Ωi \Γ for i = 1, 2.
conjugate heat transfer problem. By modifying the standard approach based on variational analysis, residuals and the generalized Green’s function, we derive accurate error estimates and use these estimates to guide adaptive mesh refinement. Modifications to the standard error analysis account for the transfer of error between components of the decomposed operator and interpolation error. Our analysis accounts for the differences between the adjoints of the full problem and the operator decomposition system. We show that the loss of order typically observed in the operator decomposition method is due to inaccuracy of transferred gradient information and we show how the boundary flux method can be used to efficiently regain the expected optimal order of convergence. REFERENCES [1] I. Babuska. The finite element method for elliptic equations with discontinuous coefficients. Computing, 5:207–213, 1970. [2] Ivo Babuska and Manil Suri. The p and hp versions of the finite element method, basic principles and properties. SIAM Review, 36:578–632, 1994. [3] Wolfgang Bangerth and Rolf Rannacher. Adaptive Finite Element Methods for Differential Equations. Birkhauser Verlag, 2003. [4] R. Becker and R. Rannacher. An optimal control approach to a posteriori error estimation in finite element methods. Acta Numerica, pages 1–102, 2001. [5] James H. Bramble and J. Thomas King. A finite element method for interface problems in domains with smooth boundaries and interfaces. Advances in Comp. Math., 67:1–19, 1998. [6] T. Cao, D.W. Kelly, and M. Ainsworth. Some useful techniques for pointwise and local error estimates of the quantities of interest in the finite element approximation. ANZIAM, 42:317–339, 2000. [7] G.F. Carey. Derivative calculation from finite element solutions. Comp. Meth. in Applied Mech. and Engr., 35:1–14, 1982. [8] V. Carey, D. Estep, and S. Tavener. A-posteriori error control of one-way coupled elliptic systems. In preparation, 2006. [9] Zhiming Chen and Jun Zou. Finite element methods and their convergence for elliptic and parabolic interface problems. Numer. Math, 79:175–202, 1998. [10] K. Eriksson, D. Estep, P. Hansbo, and C. Johnson. Introduction to adaptive methods for differential equations. Acta Numerica, pages 105–158, 1995. [11] K. Eriksson, D. Estep, P. Hansbo, and C. Johnson. Computational Differential Equations. Cambridge University Press, New York, 1996. [12] D. Estep. A posteriori error bounds and global error control for approximations of ordinary differential equations. SIAM Journal on Numerical Analysis, 32:1–48, 1995. [13] Donald Estep. A short course on duality, adjoint operators, green’s functions, and a-posteriori
CONJUGATE HEAT TRANSFER
21
error analysis. Sandia National Laboratories, Albuquerque, New Mexico, 2004. [14] Donald Estep, Michael Holst, and Mats Larson. Generalized Green’s functions and the effective domain of influence. SIAM Jour. Sci. Comp., 2004. [15] Donald Estep, Roy Williams, and Mats Larson. Estimating the error of numerical solutions of systems of reaction-diffusion equations. Memoirs of the American Mathematical Society, 146, 2000. [16] M.K. Seager G.F. Carey, S.S. Chow. Approximate boundary-flux calculations. Comp. Meth. in Applied Mech. and Engr., 50:107–120, 1985. [17] M. Giles, M. G. Larson, J. M. Levenstam, and E. Suli. Adaptive error control for finite element approximations of the lift and drag coefficients in viscous flow. Technical Report NA-97-06, Oxford University Computing Laboratory, 1997. [18] M.B. Giles. Stability analysis of numerical interface boundary conditions in fluid-structure thermal analysis. International Journal for Numerical Methods in Fluids, 25:421–436, 1997. [19] V. Ginting, D. Estep, J. Shadid, and S. Tavener. A-posteriori analysis of operator splitting for ordinary differential equations. In preparation, 2006. [20] V. Heuveline and R. Rannacher. Duality-based adaptivity in the hp-finite element method. J. Numer. Math, 0:1–18, 2003. [21] Cornelius Lanczos. Linear Differential Operators. Dover Publications, 1997. [22] Ahmed Naga and Zhimin Zhang. A posteriori error estimates based on the polynomial preserving recovery. SIAM J. Numer. Anal., 42:1780–1800, 2004. [23] J.R. Rice, P. Tsompanopoulus, and E. Vavalis. Interface relaxation methods for elliptic differential equations. Applied Numerical Mathematics, 32:219–245, 2000. [24] Barry Smith, Petter Bjorstad, and William Gropp. Domain Decomposition. Cambridge University Press, 1994. [25] M.F. Wheeler. A Galerkin procedure for estimating the flux for two-point boundary-value problems using continuous picewise-polynomial spaces. Numer. Math, 2:99–109, 1974. [26] T. Wildey. A posteriori analysis of operator decomposition in interface problems. Ph.D. Thesis, Colorado State University, 2007. [27] T. Wildey, D. Estep, and S. Tavener. A posteriori error estimation of approximate boundary fluxes. In preparation, 2006. [28] Daoqi Yang. A parallel nonoverlapping schwarz domain decomposition method for elliptic interface problems. IMA Journal on Numerical Analysis, 16:75–91, 1996. [29] Ivan Yotov. A multilevel newton-krylov interface solver for multiphysics coupling of flow in porous media. Numer. Linear Algebra Appl., 8:551–570, 2001. [30] O.C. Zienkiewicz and J.Z. Zhu. The superconvergent patch recovery and a posteriori error estimates. Intert. J. Numer. Methods Engrg., 33:1331–1364 and 1365–1382, 1992.
1. Introduction. In this paper, we consider the solution of a conjugate heat transfer problem by an operator decomposition approach. The goal is to compute a functional of the temperature of a body composed of two distinct components that share a common boundary or interface. The two components may have different conductivities and be subject to different heat sources and boundary conditions. The model consists of a system of boundary value problems for two elliptic equations coupled through boundary conditions posed on the common interface. Our interest is in situations where there is a small number of interfaces between the components. Operator decomposition provides an alternative strategy to solving the fully coupled system. We compute the temperature in each component individually subject to the conditions on the common interface that are obtained from a numerical solution of the other component. Obtaining a solution of the full system by operator decomposition therefore requires a (nominally) infinite iteration in which updated boundary conditions are passed between components. Operator decomposition is a widely used technique for solving multi-physics, multi-scale problems. The general approach is to decompose the problem into components involving simpler physics over a relatively limited range of scales, and then to seek the solution of the entire system through an iterative procedure involving solutions of the individual components. This approach is appealing because there is ∗ 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), the National Aeronautics and Space Administration (NNG04GH63G), the National Science Foundation (DMS-0107832, DGE0221595003, MSPA-CSE-0434354) and the Sandia Corporation (PO299784) † Department of Mathematics, Colorado State University, Fort Collins, CO 80523 ([email protected]). S. Taveners’s work is supported in part by the Department of Energy (DE-FG02-04ER25620) and the National Science Foundation (xxx) ‡ Department of Mathematics, Colorado State University, Fort Collins, CO 80523 ([email protected]). T. Wildey’s work is supported in part by the Department of Energy (DE-FG02-04ER25620), the National Science Foundation (DMS-0107832) and the Sandia Corporation (PO299784)
1
2
D. ESTEP, S. TAVENER, AND T. WILDEY
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. In the case of conjugate heat transfer, operator decomposition can be viewed as domain decomposition that allows for very different discretizations in each component. However, 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. In the case of the conjugate heat transfer, the operator decomposition causes a loss in the approximation order of convergence. In this paper, we perform an a posteriori error analysis of a finite-element implementation of the operator decomposition technique and obtain accurate error estimates that are used to guide an adaptive discretization strategy. Our approach is based on the standard techniques using variational analysis, residuals and the generalized Green’s function solution to an adjoint problem [3, 6, 14, 15, 12, 20], which we modify to account for several new features arising from the operator decomposition. These include the following. • Numerical errors in the solution of each component are propagated through the boundary condition applied to the other component. • Numerical errors in the solution of each the components at one step of the iterative procedure are propagated to the next step. • The adjoint of the full problem and the adjoint of the iterative procedure are distinct operators whose differences must be recognized when seeking to obtain accurate error estimates. These kinds of effects are characteristic of operator decomposition discretizations, e.g. [8, 19], and generally require extensions of the usual a posteriori analysis techniques. In addition to obtaining accurate estimates, we seek to improve the accuracy of the operator decomposition method in an efficient way. In particular, we adapt the “boundary element flux” technique developed by Wheeler [25] and Carey [16, 7] to compute normal derivatives on a boundary, and show that this can be used to improve accuracy, and in particular, restore the order of convergence that is lost due to the operator decomposition. The operator decomposition method we consider is closely related to the nonoverlapping domain decomposition methods as described in, for example, [23, 28, 24, 29]. The major concern of these authors is however the convergence of different iterative methods and the design of effective pre-conditioners or relaxation parameters, while the accuracy of the final finite element approximation is usually not addressed, despite the fact that some of the numerical results indicate a loss of accuracy. To our knowledge, this paper is the first to address the loss of accuracy. Our results would extend to this particular type of non-overlapping domain decomposition. In §2, we introduce the conjugate heat transfer problem and provide some notation. We describe the iterative operator decomposition finite element method and some modifications in §3, as well as the boundary flux method used to compute gradients on the common interface. We perform two a posteriori error analyses in §4, using first the adjoint to the fully coupled problem and then the adjoint to the iterative scheme, and present numerical results to highlight the differences between the two. We end the section with a brief discussion of adaptive mesh refinement. In §5, we carry out an analysis which identifies the transferred gradient information as being responsible for the loss of order and show that using the boundary flux method to compute the gradient information restores the order of convergence. Our conclusions
3
CONJUGATE HEAT TRANSFER
are presented in §6. 2. The model for conjugate heat transfer. Let Ω1 and Ω2 be polygonal domains in R2 or R3 with boundaries ∂Ω1 and ∂Ω2 intersecting along an interface Γ = ∂Ω1 ∩ ∂Ω2 . We consider a system of second order linear elliptic problems, where the components are coupled through boundary conditions imposed on Γ, L1 u1 = f1 , u (1 = 0, u =u , 1 2 A1 ∂n u1 = A2 ∂n u2 , L2 u2 = f2 , u2 = 0,
x ∈ Ω1 , x ∈ ∂Ω1 \Γ, x ∈ Γ,
(2.1)
x ∈ Ω2 , x ∈ ∂Ω2 \Γ,
where for i = 1, 2, Li ui = −∇ · (Ai ∇ui ) + ci ui , Ai ≥ Ai,0 > 0, ci , fi are sufficiently smooth functions and ∂n denotes the unit outward normal derivative to ∂Ω1 . The results of this paper extend easily to general elliptic operators and general Dirichlet, Neumann, and Robin boundary conditions on the boundaries in the complement of the interface. We let L2 (Ωi ) denote the space of square integrable functions on Ωi with inner product (·, ·)Ωi and norm k · kΩi , but use (·, ·) = (·, ·)Ωi when the domain is clear. We use H s (Ωi ) to denote the Sobolev space with real index s associated with the norm k · kΩi ,s and seminorm | · |Ωi ,s . We also use the subspaces H01 (Ωi ) = 1 v ∈ H (Ωi ), v = 0 on ∂Ωi \Γ . The weak formulation of (2.1) seeks ui ∈ H01 (Ωi ) such that u1 = u2 on Γ and 2 X i=1
ai (ui , vi ) =
2 X
(fi , vi ),
(2.2)
i=1
R for all vi ∈ H01 (Ωi ) with ai (ui , v) = Ωi (Ai ∇ui · ∇v + ci ui v) dx for i = 1, 2. Assuming that each ai (·, ·) is coercive, (2.2) admits a unique weak solution in H01 (Ωi ) [1, 5, 9]. 3. An iterative operator decomposition method. To compute a numerical (0) solution of (2.1), we consider iterative operator decomposition methods. Let u2 be an initial guess for the Dirichlet data along the interface. We solve (k) x ∈ Ω1 , L1 u1 = f1 , (k) (3.1) u1 = 0, x ∈ ∂Ω1 \Γ, (k) (k−1) u1 = u2 , x ∈ Γ, followed by
(k) L2 u2 = f2 , (k) u2 = 0, (k) (k) A1 ∂n u1 = A2 ∂n u2 ,
x ∈ Ω2 , x ∈ ∂Ω2 \Γ, x ∈ Γ.
(3.2)
(k)
We iterate until (hopefully) a convergence criteria is met, e.g., ku2 smaller than a prescribed tolerance.
(k−1)
− u2
kΓ is
4
D. ESTEP, S. TAVENER, AND T. WILDEY
3.1. Finite element discretization. We let Ti,h be a triangulation of Ωi into elements K where the length of the longest edge is hK and h = maxK∈Ti,h hK . We assume that each triangulation is locally quasi-uniform and Ωi = ∪K∈Ti,h K. However, the triangulations on either side of Γ are not assumed to be aligned. We use the piecewise polynomial spaces S1 = v continuous on Ω1 , v ∈ P 1 (K) for all K ∈ T1,h , S2 = v continuous on Ω2 , v ∈ P 1 (K) for all K ∈ T2,h , and the associated spaces
S1,0 = {v ∈ S1 | v = 0, x ∈ ∂Ω1 } , S2,0 = {v ∈ S2 | v = 0, x ∈ ∂Ω2 \Γ} , where P 1 (K) denotes the space of linear polynomials on an element K. We let πi be a projection into Si as well as the projection into Si along the interface Γ. (k) For the finite element approximation, we compute U1 ∈ S1 satisfying (k) a1 (U1 , v) = (f1 , v)Ω1 , for all v ∈ S1,0 , (k) (3.3) x ∈ ∂Ω1 \Γ, U1 = 0, (k) (k−1) U1 = π1 U2 , x ∈ Γ, (k)
followed by U2 ∈ S2 such that ( (k) (k) a2 (U2 , v) = (f2 , v)Ω2 − (A1 ∂n U1 , v)Γ , (k) U2 = 0,
for all v ∈ S2,0 , x ∈ ∂Ω2 \Γ.
(3.4)
3.2. Relaxed iterations. Unfortunately, the simple iterative (3.3)-(3.4) may not converge. In particular, the convergence depends on the values of A1 and A2 along the interface and the geometry of each region [18, 23, 28]. As an alternative, we consider two “relaxed” iteration schemes. (1) Relaxed Dirichlet values We choose α ∈ [0, 1] and update the Dirichlet interface values with (k)
U1
(k−1)
= αU1
(k−1)
+ (1 − α)π1 U2
.
(3.5)
Optimal values of α can be found in [23, 28], but α = 1/2 works well in most situations. (2) Relaxed Neumann values (k−1) (k) We chose β ∈ [0, 1], set Nβ = βA2 ∂n U2 + (1 − β)A1 ∂n U1 , and seek (k) U2 ∈ S2 such that (k)
a2 (U2 , v) = (f2 , v)Ω2 − (Nβ , v)Γ , for all v ∈ S2,0 ,
(3.6)
A proper choice of β reduces the number of iterations. 3.3. Flux correction. It turns out that the operator decomposition results in a loss of order arising from the relatively low order of accuracy in the derivative of the finite element solution, which is O(h) when computed with a standard O(h2 ) method. We show below that this reduces the order of the overall approximation to first order.
CONJUGATE HEAT TRANSFER
5
To mitigate this effect, we use a post-processing technique developed by Wheeler [25] and Carey [16, 7] to compute more accurate boundary flux information. We define the set of elements in T1,h that intersect the boundary by Γ T1,h = {K ∈ T1,h | K ∩ ∂Ω 6= ∅},
and the corresponding space Γ Wh = {v ∈ P 1 (K) with K ∈ T1,h , v(ηi ) = 0 if ηi ∈ / ∂Ω},
where {ηi } denotes the nodes of element K, so the degrees of freedom correspond to the nodes on the boundary. We seek σ (k) ∈ Wh satisfying (k)
−(σ (k) , v)Γ = (f1 , v)Ω1 − a1 (U1 , v), for all v ∈ Wh ,
(3.7)
(k)
where U1 solves (3.3). Green’s identity implies that σ (k) gives an approximation to the normal flux on the boundary which is relatively inexpensive to compute. In general, the accuracy of the flux approximation depends on the regularity of an associated Green’s function [17, 27], but in many cases this post-processed flux is O(h2 ) rather than the standard O(h) for the normal flux of a piecewise linear finite element approximation. It turns out that a fortunate cancellation of errors makes the accuracy of this recovered flux unimportant for our purposes. We stress that σ (k) is not assumed to be continuous. In fact, if the domain has a corner on the boundary, the normal derivative is, in general, discontinuous due to the jump in the normal vector. When Dirichlet conditions are given on each boundary segment we implement the method described in [7] and allow two degrees of freedom to account for the discontinuity. For the application to more general boundary conditions, such as a Neumann or Robin condition on one segment, or an interface condition on both, we refer the reader to [26]. Now, as an alternative to (3.4), we may solve ( (k) a2 (U2 , v) = (f2 , v)Ω2 − (σ (k) , v)Γ , for all v ∈ S2,0 , (3.8) (k) U2 = 0, x ∈ ∂Ω2 \Γ. Other possible approaches to mitigating the loss of order include approximating the boundary flux using a gradient recovery techniques such as the ZienkiewiczZhu (ZZ) patch recovery technique [30] or the polynomial preserving recovery (PPR) method [22], or using higher order polynomials near the interface to improve the accuracy of the finite element flux [2, 20]. We had mixed computational success using these other approaches. 4. A posteriori error analysis. To estimate the error of the operator decomposition finite element approximation, we apply a posteriori techniques based on variational analysis and the adjoint problem. In this case, however, the adjoint for the fully coupled original problem differs significantly from the adjoint problem associated with the discretization of the decomposed system. The motivation to use operator decomposition suggests that the adjoint of the full problem will be unavailable in practice and therefore the decomposed adjoint will be used to compute error estimates. The difference between the adjoint problems can be decreased by increasing the computational work used to solve the decomposed adjoint, however in practice, we must consider the issue that the two adjoints lead to significantly different error representations.
6
D. ESTEP, S. TAVENER, AND T. WILDEY
In order to understand the differences, we derive a posteriori error estimates using both adjoints. Both analyses begin in the same way. We wish to estimate the difference between the exact solution to the full problem and the numerical approximation to the solution of the iterative procedure, i.e. E (k) = u − U (k) . We decompose the error, E (k) = u − U (k) = u − u(k) + u(k) − U (k) = c(k) + e(k) .
The first component c(k) measures the difference between the exact solution to the full problem and the exact solution to the iterative problem. If the iterative method converges, then c(k) → 0. The second component e(k) measures the numerical error in solving the iterative problem. It is helpful to remark that in the conventional approach, the adjoint is defined in such a way as to make the formal bilinear identity hold [21, 13], i.e., to make the boundary terms arising from integration by parts equal to zero. In the analysis below, we find it convenient to consider the values passed between components as quantities of interest defined on the common boundary and we have to alter the definition of the adjoint problem accordingly. 4.1. The adjoint to the fully coupled problem. The adjoint boundary value problem for the quantity of interest (ψ, u) = (ψ1 , u1 )+(ψ2 , u2 ) for the coupled problem (2.1) is ∗ x ∈ Ω1 , L1 φ1 = ψ1 , φ x ∈ ∂Ω1 \Γ, (1 = 0, φ =φ , 1 2 (4.1) x ∈ Γ, A ∂ φ 1 n 1 = A2 ∂n φ2 , L∗2 φ2 = ψ2 , x ∈ Ω2 , φ2 = 0, x ∈ ∂Ω2 \Γ, where L∗i φi = −∇ · (Ai ∇φi ) + ci φi . We use a∗i (·, ·) to represent the weak form of the adjoint operator. We solve (4.1) numerically by using an iterative operator decomposition approach as for the forward problem. The iterations are completely independent of the forward iterations. We can derive an error representation formula for the basic scheme (3.1)-(3.2), a weighted relaxation technique (3.5) or (3.6), or when using the post-processed flux as (k) in (3.8). In the discussion below, we use θh to denote the numerical flux passed at the k th iteration from Ω1 to Ω2 . To begin, we multiply the system (4.1) by (ψ1 , ψ2 )T and apply the divergence theorem, noting that u1 = u2 and A1 ∂n φ1 = A2 ∂n φ2 along Γ, to obtain (ψ1 , e1 ) + (ψ2 , e2 ) = a1 (e1 , φ1 ) + a2 (e2 , φ2 ) (k)
(k)
+ (U1 , A1 ∂n φ1 )Γ − (U2 , A2 ∂n φ2 )Γ . Observe that the test space S1,0 consists of functions that are zero along the interface, while in general, φ1 is not zero along Γ. This means that the projection of φ1 into S1,0 cannot be the interpolant. We define a new projection π10 : H 2 → S1,0 such that for any node xi ( π1 φ(xi ), xi ∈ / Γ, π10 φ(xi ) = (4.2) 0, xi ∈ Γ.
7
CONJUGATE HEAT TRANSFER
We also observe that (k)
a2 (e2 , π2 φ2 ) = −(A1 ∂n u1 , π2 φ2 )Γ + (θh , π2 φ2 )Γ , and a1 (u1 , φ1 − π10 φ1 ) = (f1 , φ1 − π10 φ1 ) + (A1 ∂n u1 , φ1 )Γ , a2 (u2 , φ2 − π2 φ2 ) = (f2 , φ2 − π2 φ2 ) − (A1 ∂n u1 , φ2 − π2 φ2 )Γ , since the adjoint solutions are not zero along Γ. Using the projection (4.2) in the Galerkin orthogonality relation and the equalities above, we have (k)
(ψ1 , e1 ) + (ψ2 , e2 ) = (f1 , φ1 − π10 φ1 ) − a1 (U1 , φ1 − π10 φ1 ) (k)
+ (f2 , φ2 − π2 φ2 ) − a2 (U2 , φ2 − π2 φ2 ) (k)
+ (U1
(k)
(k)
− U2 , A1 ∂n φ1 )Γ + (θh , π2 φ2 )Γ .
Next, we define π∂ φ1 = π1 φ1 − π10 φ1 which is nonzero only near the interface due to the definition of π10 φ1 . Substituting π10 φ1 = π1 φ1 − π∂ φ1 gives (k)
(ψ1 , e1 ) + (ψ2 , e2 ) = (f1 , φ1 − π1 φ1 ) − a1 (U1 , φ1 − π1 φ1 ) (k)
+ (f2 , φ2 − π2 φ2 ) − a2 (U2 , φ2 − π2 φ2 ) (k)
+ (f1 , π∂ φ1 ) − a1 (U1 , π∂ φ1 ) (k)
+ (U1
(k)
(k)
− U2 , A1 ∂n φ1 )Γ + (θh , π2 φ2 )Γ
Finally (3.7) implies that the recovered boundary flux, σ (k) satisfies (k)
−(σ (k) , π∂ φ1 )Γ = (f1 , π∂ φ1 )Ω1 − a1 (U1 , π∂ φ1 ), while π∂ φ1 = π1 φ1 along Γ. We conclude (k) (k) Theorem 4.1. The errors e1 = u1 − U1 and e2 = u2 − U2 satisfy (k)
(ψ1 , e1 ) + (ψ2 , e2 ) = (f1 , φ1 − π1 φ1 ) − a1 (U1 , φ1 − π1 φ1 ) (k)
+ (f2 , φ2 − π2 φ2 ) − a2 (U2 , φ2 − π2 φ2 ) (k)
+ (U1
(k)
− U2 , A1 ∂n φ1 )Γ
(k)
+ (θh , π2 φ2 )Γ − (σ (k) , π1 φ1 )Γ
(4.3)
The error has been decomposed into two typical discretization components, an iterative component, and a component reflecting the error arising from the transfer of derivative information. The choice of iterative method does not affect the discretization component of the error, but it does influence both the iterative and the transfer components along Γ. We illustrate this for a few common iterative schemes: (k) (k−1) • Suppose U1 = π1 U2 , then (k)
(U1
(k)
− U2 , A1 ∂n φ1 )Γ (k−1)
= (π1 U2
(k−1)
− U2
(k)
, A1 ∂n φ1 )Γ + (U2
(k−1)
− U2
which represents a projection error and an iteration error.
, A1 ∂n φ1 )Γ ,
8
D. ESTEP, S. TAVENER, AND T. WILDEY (k)
• Suppose U1 (k)
(U1
(k−1)
= αU1
(k−1)
+ (1 − α)π1 U2
(k)
− U2 , A1 ∂n φ1 )Γ =
k−2 X
, then (k−i−1)
αi−1 (π1 (U1
(k−i)
− U2
), A1 ∂n φ1 )Γ
i=1
(k−1)
+ (U2
(k)
− U2 , A1 ∂n φ1 )Γ
(k−1)
+ (π1 U2
(k−1)
− U2
, A1 ∂n φ1 )Γ ,
which represents a series of iteration errors and a projection error. Notice that since α < 1, the effect of the iteration error from previous iterations diminishes due to the increasing power on α. (k) (k) • Suppose we set θh = A1 ∂n U1 , then (k)
(k)
(θh , π2 φ2 )Γ −(σ (k) , π1 φ1 )Γ = (A1 ∂n U1 −σ (k) , π2 φ2 )Γ +(σ (k) , π2 φ2 −π1 φ1 )Γ , which represents a transfer error and a projection error. (k) • Suppose we set θh = σ (k) , then (k)
(θh , π2 φ2 )Γ − (σ (k) , π1 φ1 )Γ = (σ (k) , π2 φ2 − π1 φ1 )Γ , which represents only a projection error with no transfer error. 4.2. The adjoint to the iterative scheme. Next, we derive an error representation using the natural adjoint for the iterative system (3.1)-(3.2), where we set (k) (k−1) (k) (k) U1 = π1 U2 and A2 ∂n U2 = A1 ∂n U1 . This adjoint avoids the formulation of a globally coupled adjoint, and reads (k) ∗ (k) L1 φ1 = ψ1 , x ∈ Ω1 , (k) φ1 = 0, x ∈ ∂Ω1 \Γ, (k) (k) φ1 = φ2 , x ∈ Γ, (k) ∗ (k) L2 φ2 = ψ2 , (k) φ2 = 0, (k) (k+1) A2 ∂n φ2 = A1 ∂n φ1 ,
x ∈ Ω2 , x ∈ ∂Ω2 \Γ, x ∈ Γ,
(N +1)
(4.4)
(4.5)
for k = N, . . . , 1 with A1 ∂n φ1 = 0. Note that the adjoint system is defined “backwards”. In the reasonable case in which we seek information from the final (1) (1) (N −1) iterate, the data for the adjoint problem would be ψ1 = ψ2 = . . . = ψ1 = (N −1) (N ) (N ) ψ2 = 0, with ψ1 and ψ2 chosen appropriately. We derive the error representation formula for this system by observing N X
k=1
N X (k) (k) (k) (k) (k) (k) (k) (k) (L∗1 φ1 , e1 ) + (L∗2 φ2 , e2 ) , (ψ1 , e1 ) + (ψ2 , e2 ) = k=1
and apply the steps used to derive (4.3). This provides the following theorem.
9
CONJUGATE HEAT TRANSFER (1)
(1)
(N )
(N )
Theorem 4.2. The errors e1 , e2 , . . . , e1 , e2 N X
k=1
satisfy
(k) (k) (k) (k) (ψ1 , e1 ) + (ψ2 , e2 )
=
N X
(k)
(k)
(k)
(k)
(k)
(f1 , φ1 − π1 φ1 ) − a1 (U1 , φ1 − π1 φ1 )
k=1
(k)
(k)
(k)
(k)
(k)
+ (f2 , φ2 − π2 φ2 ) − a2 (U2 , φ2 − π2 φ2 ) +
(k) (U1
−
(k−1) (k) U2 , A1 ∂n φ1 )Γ
+
(k) (k) (A1 ∂n U1 , π1 φ1 )Γ
− (σ
(k)
(k) , π2 φ2 )Γ
(4.6) It is possible to define adjoints for more complicated iterative methods, but it (k) (k−1) may be difficult. For example, relaxation techniques couple U1 to both U2 and (k−1) U1 which affects how the adjoint variables are coupled. In addition, using the post-processed numerical flux given in section 3.2 requires defining the adjoint of the post-processing procedure. 4.3. Numerical results. We illustrate the a posteriori estimates in §4.1 and §4.2. Example 4.1. We triangulate the domains Ω1 = [−0.25, 0.25] × [−0.25, 0.25] and Ω2 = ([−1, 1] × [−1, 1]) \Ω1 independently, see Fig. 4.1, and consider the elliptic
Fig. 4.1. Triangulations of Ω1 and Ω2 that do not match along the interface.
interface problem, −∇ · (∇u1 ) = f1 , ( u1 = u2 , ∂n u 1 = ∂n u 2 , −∇ · (∇u2 ) = f2 , u2 = 0,
x ∈ Ω1 , x ∈ Γ,
(4.7)
x ∈ Ω2 , x ∈ ∂Ω2 \Γ,
where f1 (x, y) and f2 (x, y) are chosen so the true solutions are u2 = sin(4πx) sin(4πy) (k) (k−1) and u1 = 3u2 . We solve this using (3.3)-(3.4) with U1 = π1 U2 and iterate until (k) (k−1) −5 kU1 − U1 kΓ < 10 . We consider several quantities of interest
10
D. ESTEP, S. TAVENER, AND T. WILDEY
1. 2. 3. 4.
the global average value, the average value over Ω1 , the average value over Ω2 , the value of u2 at the point (0, 0) computed using the approximate delta 2 2 function δˆ0 (x, y) = 400 π exp(−400x − 400y ). First we use an independent iterative scheme to solve the adjoint of the fully (k) coupled problem using the relaxation parameter α = 1/2 and iterating until kΦ1 − (k−1) (k) Φ1 kΓ < 10−6 , where Φ1 represents the approximation of the adjoint solution at th the k iteration. We display the error estimates for each of the quantities of interest in Table 4.1. Next we solve the adjoint to the operator decomposition method and give these results in Table 4.2. Note that the number of iterations, as well as the values of the relaxation parameters, are determined by the number of iterations used in the solution of the forward operator decomposition method.
1. 2. 3. 4.
ψ1 4 4 0 δˆ0
ψ2 4/15 0 4/15 0
True Error −0.20 −0.15 −0.046 −0.15
F.C. Adjoint −0.20 −0.15 −0.046 −0.15
Effect. Ratio 0.99 0.99 1.00 0.99
Table 4.1 Error estimates and effectivity ratios using the adjoint for the fully coupled system.
1. 2. 3. 4.
ψ1 4 4 0 δˆ0
ψ2 4/15 0 4/15 0
True Error −0.20 −0.15 −0.046 −0.15
O.D. Adjoint −0.20 −0.15 −0.046 −0.15
Effect. Ratio 0.99 0.99 1.00 0.99
Table 4.2 Error estimates and effectivity ratios using the adjoint for the operator decomposition method.
Example 4.2. Next we consider neighboring domains Ω1 = [0, 1]×[0, 1] and Ω2 = ([1, 2] × [0, 1]) with independent triangulations, and the elliptic interface problem −∇ · (A1 ∇u1 ) = f1 , u (1 = 0, u =u , 1 2 A ∂ u 1 n 1 = A2 ∂n u2 , −∇ · (A2 ∇u2 ) = f2 , u2 = 0,
x ∈ Ω1 , x ∈ ∂Ω1 \Γ, x ∈ Γ,
(4.8)
x ∈ Ω2 , x ∈ ∂Ω2 \Γ,
where A1 = 1, A2 = 3, and f1 (x, y) and f2 (x, y) are chosen so the true solutions are u2 = sin(2πx) sin(2πy) and u1 = 3u2 . The quantity of interest is the average value over Ω1 ∪ Ω2 , so ψ1 = ψ2 = 1. This experiment demonstrates that the error representation formulas (4.3) and (4.6) may give different estimates if relatively few iterations are used in the operator
11
CONJUGATE HEAT TRANSFER 1.2 1.6
Adjoint for Fully Coupled Problem Adjoint for Operator Decomposition Method
Size of adjoints
Effectivity Ratios
1.1
1
0.9
1.2
0.8
0.4 0.8 2
4
6
8 10 Iterations
12
14
0
16
0
5
10
15
Iteration
Fig. 4.2. On the left, a comparison of the effectivity ratios using the two adjoints for a given number of iterations used for the operator decomposition method for the forward problem. On the (k) (k) right, we plot kφ1 k+kφ2 k to show the decay of influence of errors that occur in previous iterations.
decomposition method for the forward problem. We set a fixed number of iterations in the operator decomposition method for (4.8) and use (3.3) and (3.4) to approximate the solution. We solve the fully coupled adjoint iteratively using sufficient iterations for convergence. However, the number of iterations for the adjoint of the iterative system is determined by the number of iterations for the forward problem. In Fig. 4.2, we plot the number of iterations versus the effectivity ratios for the two estimates. We observe that the adjoint for the operator decomposition method does not produce accurate estimates until a sufficient number of iterations for the forward problem have been carried out. This implies that the iterative component of the error, c(k) = u − u(k) , dominates the estimate. If the iterative method converges, we expect that the effect of the initial guess and the first few iterations would diminish. This is reflected in Fig. 4.2 where we (k) (k) plot the L2 norm of the iterative adjoint φ1 and φ2 over Ω1 and Ω2 respectively at each iteration. We observe that the norm decays rapidly as k decreases, which indicates that the norm of the fixed point operator, and hence the norm of the adjoint fixed point operator, is small. Thus, the influence of errors in previous iterations are rapidly damped and have little influence on the current error. This implies that we can compute accurate estimates using the truncated series (N )
(ψ1
(N )
(N )
(N )
, e1 ) + (ψ2 , e2 ) N X (k) (k) (k) (k) (k) (f1 , φ1 − π1 φ1 ) − a1 (U1 , φ1 − π1 φ1 ) = k=M
(k)
(k)
(k)
(k)
(k)
+ (f2 , φ2 − π2 φ2 ) − a2 (U2 , φ2 − π2 φ2 ) (k)
+ (U1
(k−1)
− U2
(k)
(k)
(k)
(k)
, A1 ∂n φ1 )Γ + (A1 ∂n U1 , π1 φ1 )Γ − (σ (k) , π2 φ2 )Γ
(4.9) where 1 ≤ M ≤ N , in place of the full series (4.6). This significantly reduces the computational cost associated with computing a series of adjoint problems. In Fig. 4.3 we plot the effectivity ratio for the global average value where we truncate the error representation and compute only the last N terms.
12
D. ESTEP, S. TAVENER, AND T. WILDEY
Effectivity Ratio
1.15 1.05 0.95 0.85 0.75 0.65 0
2
4
6 8 10 Number of Terms
12
14
Fig. 4.3. A comparison of the effectivity ratios computed using a truncated error representation.
4.4. Adaptive refinement. We use the a posteriori error estimate as the basis for adaptivity by employing the standard “optimization framework” after writing the estimate as a sum of element contributions and introducing norms [10, 11, 4, 3]. An element K is marked for refinement when the local error indicator is larger than a local tolerance, usually the global tolerance divided by the current number of elements. Example 4.3. We demonstrate the adaptive procedure on the model problem (4.7) with the quantity of interest equal to the value of u2 at the point (1.75, 0.25), which we approximate by choosing ψ1 = 0 and ψ2 = 400/π exp(−400(x − 1.75)2 − 400(y − 0.25)2 ). Fig. 4.4 shows the results after 3 refinement steps. We solve (3.3)(3.4) and observe that refinement occurs within Ω1 along Γ which reflects the impact of numerical errors in the normal derivative on the quantity of interest.
Fig. 4.4. Adaptive mesh for the quantity of interest equal to the value of u2 at the point (1.75, 0.25).
5. An analysis of the loss of order. In practice, the iterative technique (3.3) and (3.4) is occasionally observed to result in O(h) convergence rather than the O(h2 ) convergence that is obtained when solving the fully coupled problem. This loss of order results from passing the normal derivative of the finite element approximation, which is only O(h). Passing the recovered boundary flux restores the full order of convergence. We use the adjoint for the fully coupled problem to derive a posteriori (k) (k−1) error bounds for the iterative approximations in the case where U1 = π1 U2 . The case where the Dirichlet values are updated using a relaxation technique can be analyzed using the same approach and gives identical results. Numerical examples are provided at the end of the section.
CONJUGATE HEAT TRANSFER
13
5.1. L2 error bounds. Let u1 ∈ H 1+α1 (Ω1 ) and u2 ∈ H 1+α2 (Ω2 ) be the solu(k) (k) tions to (2.1) with 0 ≤ α1 , α2 ≤ 1, and U1 and U2 be the solutions of (3.3) and th 1+α1 (Ω1 ) and φ2 ∈ H 1+α2 (Ω2 ) and (3.4) respectively at the k iteration. Let φ1 ∈ H pose the adjoint problem (4.1) with ψ1 = e1 /ke1 kΩ1 and ψ2 = e2 /ke2 kΩ2 . Starting with (4.3), integration by parts over each element K gives ke1 kΩ1 + ke2 kΩ2 = I1 + I2 + I3 + I4 , where I1 =
X
K∈T1,h
1 (k) (k) (f1 − L1 U1 , φ1 − π1 φ1 )K + ([A1 ∂n U1 ], φ1 − π1 φ1 )∂K 2 X
+
K∈T2,h
1 (k) (k) (f2 − L2 U2 , φ2 − π2 φ2 )K + ([A2 ∂n U2 ], φ2 − π2 φ2 )∂K 2
(k)
(k)
I2 =(A2 ∂n U2 , φ2 − π2 φ2 )Γ − (A1 ∂n U1 , φ1 − π1 φ1 )Γ , (k)
(k)
− U2 )Γ ,
I3 =(A1 ∂n φ1 , U1 (k)
I4 =(A1 ∂n U1 , π2 φ2 )Γ − (σ (k) , π1 φ1 )Γ , where [Ai ∂n Ui ] denoting the jump in the normal derivative across an element edge. The first term I1 is a standard a posteriori error bound for elliptic problems and is not affected by non-matching triangulations along the interface or by transfer error. The second term I2 is similar to the jump terms along element edges in I1 and is the expected jump term when the triangulations align along the interface. The third term I3 represents the jump in the Dirichlet values across the interface. Finally, the fourth term I4 represents the difference between the flux passed from Ω1 to Ω2 and the flux obtained via the boundary-flux correction technique. We first construct Lemmas 5.1 to 5.4 below to bound I1 to I4 individually. In each of these Lemmas we first provide the general bound when the triangulations do not match across the boundary and then show the simplification that arises for matching triangulations. We then combine these four lemmas into Theorems 5.5 and 5.7 which give error bounds for the basic iteration (3.3) and (3.4), and when using flux correction (3.8) respectively. These two theorems describe the general result when the triangulations do not match across the boundary while the simplification given matching triangulations is provided as a Corollary. Lemma 5.1. (Bound on I1 ) I1 ≤
X
K∈T1,h
+
X
! 1 Ch1+α |φ1 |K,1+α1 K · 1 Ch1+α |φ1 |K,1+α1 K ! (k) 2 kf2 − L2 U2 kK Ch1+α |φ2 |K,1+α2 K . · −1/2 (k) 2 1 Ch1+α |φ2 |K,1+α2 [A2 ∂n U2 ]k∂K K 2 khK (k)
kf1 − L1 U1 kK −1/2 (k) 1 [A1 ∂n U1 ]k∂K 2 khK
K∈T2,h
Proof. Apply standard arguments using the Cauchy-Schwarz inequality and a trace inequality.
14
D. ESTEP, S. TAVENER, AND T. WILDEY
If u1 ∈ H 2 (Ω1 ), we can also provide a bound for the jump terms. Ordinarily, we would add and subtract (A1 ∂n u1 , φ1 − π1 φ1 )∂K to estimate the jump term. However, (k) this causes an iteration error (A1 ∂n u1 − A1 ∂n u1 , φ1 − π1 φ1 )Γ , which may be quite (k) complicated. Therefore, we let u ˜1 solve
(k)
(k) u1 , v) = (f1 , v)Ω1 , a1 (˜ (k) u ˜1 = 0, (k) (k−1) u ˜1 = π1 U2 ,
for all v ∈ H01 (Ω1 ), x ∈ ∂Ω\Γ, x ∈ Γ.
(5.1)
To be precise, u ˜1 is the exact weak solution over Ω1 given the data from the previous iteration. Standard finite element theory can be applied to estimate (k)
(k)
k˜ u1 − U1 kΩ1 ,1 ≤ Ch1 kf1 kΩ1 . (k)
We add and subtract (A1 ∂n u ˜1 , φ1 − π1 φ1 )∂K and use the Cauchy-Schwarz and trace inequalities to obtain, (k)
1/2
k[A1 ∂n U1 ]k∂K ≤ ChK kf1 kΩ1 . A similar argument holds if u2 ∈ H 2 (Ω2 ). Lemma 5.2. (Bound on I2 ) If the triangulations T1,h and T2,h do not match along the interface, then (k) −1/2 2 |φ2 |1+α2 + kA1 ∂n U1 kΓ · (kπ1 φ1 − π2 φ1 kΓ ) I2 ≤ kh2 [A∂n U (k) ]kΓ · Ch1+α 2
where
kπ1 φ1 − π2 φ1 kΓ ≤ kπ1 φ1 − φ1 kΓ + kφ1 − π2 φ1 kΓ 1/2+α1
≤ C1 h1
1/2+α2
|φ1 |Ω1 ,1+α1 + C2 h2
|φ2 |Ω2 ,1+α2 .
If the triangulations T1,h and T2,h match along the interface, then −1/2 2 |φ2 |1+α2 . I2 ≤ kh2 [A∂n U (k) ]kΓ · Ch1+α 2 (k)
where [A∂n U (k) ] = A2 ∂n U2
(k)
− A1 ∂n U1
and h2 = maxK∈T2,h hK . (k)
Proof. We add and subtract (A1 ∂n U1 , φ2 − π2 φ2 )Γ and use φ1 = φ2 on Γ to write (k)
I2 = (A2 ∂n U2
(k)
(k)
− A1 ∂n U1 , φ2 − π2 φ2 )Γ + (A1 ∂n U1 , π1 φ1 − π2 φ1 )Γ .
Observe π1 φ1 = π2 φ2 if the triangulations match along the interface. We apply the Cauchy-Schwarz and trace inequalities to complete the proof. (k)
If u2 ∈ H 2 (Ω2 ), we may define u ˜2
(k)
similarly to u ˜1 1/2
k[A∂n U ]kΓ ≤ Ch2 kf2 kΩ1 .
in (5.1) and obtain
15
CONJUGATE HEAT TRANSFER
Lemma 5.3. (Bound on I3 ) If the triangulations T1,h and T2,h do not match along the interface, then (k−1) (k) I3 ≤ (kA1 ∂n φ1 kΓ ) · kU2 − U2 kΓ (k−1) (k−1) + (kA1 ∂n φ1 kΓ ) · kπ1 U2 − U2 kΓ where
(k−1)
kπ1 U2
(k−1)
− U2
1/2+α1
kΓ ≤ C1 h1
1/2+α2
|ˆ u|Ω1 ,1+α1 + C2 h2
|ˆ u|Ω2 ,1+α2 .
If the triangulations T1,h and T2,h match along the interface, then (k−1) (k) I3 ≤ (kA1 ∂n φ1 kΓ ) · kU2 − U2 kΓ . (k)
Proof. First, observe that U1 (k−1) tracting (A1 ∂n φ1 , U2 )Γ to get (k−1)
I3 = (A1 ∂n φ1 , π1 U2 (k−1)
(k−1)
= π1 U2 (k−1)
− U2
and rewrite I3 by adding and sub(k−1)
)Γ + (A1 ∂n φ1 , U2
(k)
− U2 )Γ .
(k−1)
To bound kπ1 U2 − U2 kΓ we introduce a new function u ˆ such that u ˆ ∈ (k−1) (k−1) H 1+α1 (Ω1 ) and u ˆ ∈ H 1+α2 (Ω2 ), as well as π1 u ˆ = π1 U2 and π2 u ˆ = U2 . This yields (k−1)
kπ1 U2
(k−1)
− U2
(k−1)
1/2+α1
kΓ ≤ C1 h1
1/2+α2
|ˆ u|Ω1 ,1+α1 + C2 h2
|ˆ u|Ω2 ,1+α2 .
(k−1)
Observe that π1 U2 = U2 if the triangulations match along the interface. We apply the Cauchy-Schwarz inequality to complete the proof. Lemma 5.4. (Bound on I4 ) If the triangulations T1,h and T2,h do not match along the interface, then (k) (k) (k) I4 ≤ kA1 ∂n U1 − A1 ∂n u ˜1 kΓ + kA1 ∂n u ˜1 − σ (k) kΓ · (kπ2 φ2 kΓ ) + kσ (k) kΓ · (kπ2 φ2 − π1 φ1 kΓ )
where
1/2+α1
kπ2 φ2 − π1 φ1 kΓ ≤ C1 h1
1/2+α2
|φ1 |Ω1 ,1+α1 + C2 h2
|φ2 |Ω2 ,1+α2 .
If the triangulations T1,h and T2,h match along the interface, then (k) (k) (k) I4 ≤ kA1 ∂n U1 − A1 ∂n u ˜1 kΓ + kA1 ∂n u ˜1 − σ (k) kΓ · (kπ2 φ2 kΓ ) . (k)
where u ˜1
is defined by (5.1).
Proof. We add and subtract (σ (k) , π2 φ2 )Γ and use φ1 = φ2 along Γ to get (k)
I4 = (A1 ∂n U1
− σ (k) , π2 φ2 )Γ + (σ (k) , π2 φ2 − π1 φ2 )Γ .
16
D. ESTEP, S. TAVENER, AND T. WILDEY
Observe π1 φ1 = π2 φ2 if the triangulations match along the interface. We add and (k) subtract (A1 ∂n u ˜1 , π2 φ2 )Γ to the first term and apply the Cauchy-Schwarz inequality to complete the proof.
In practice, the error in the normal derivative is typically the same accuracy as 1 the H 1 error, namely O(hα 1 ). However, an application of the trace theorem only α1 /2 proves O(h1 ) accuracy. This is not an important issue, however, since we intend to use the fact that this term is less accurate than the others, and therefore pollutes the L2 error. We assume the error in the normal derivative can be bounded, (k)
kA1 ∂n U1
(k)
−u ˜1 kΓ ≤ Chβ1 kf1 kΩ1 ,
for α1 /2 ≤ β ≤ α1 . The error in the recovered boundary flux can be bounded (k)
1 kf1 kΩ1 , kA1 ∂n u ˜1 − σ (k) kΓ ≤ CS1 h1+α 1
where S1 is a stability factor defined by an associated Green’s function [17, 27].
Theorem 5.5. Assume the triangulations T1,h and T2,h do not match along (k) (k) the interface Γ. If U1 and U2 solve (3.3) and (3.4) respectively, then the errors (k) (k) e1 = u1 − U1 and e2 = u2 − U2 satisfy
ke1 kΩ1 + ke2 kΩ2 ≤
X
K∈T1,h
+
X
K∈T2,h
! 1 Ch1+α |φ1 |K,1+α1 K · 1 Ch1+α |φ1 |K,1+α1 K ! (k) 2 kf2 − L2 U2 kK Ch1+α |φ2 |K,1+α2 K · −1/2 (k) 2 1 Ch1+α |φ2 |K,1+α2 [A2 ∂n U2 ]k∂K K 2 khK −1/2 2 |φ2 |1+α2 + kh2 [A∂n U (k) ]kΓ · Ch1+α 2 (k−1) (k) + (kA1 ∂n φ1 kΓ ) · kU2 − U2 kΓ 1 + Chβ1 kf1 kΩ1 + CS1 h1+α kf1 kΩ1 · (kπ2 φ2 kΓ ) 1 (k)
kf1 − L1 U1 kK −1/2 (k) 1 [A1 ∂n U1 ]k∂K 2 khK
+C1 h1
1/2+α1
u|Ω1 ,1+α1 ) ((S2 + S4 )|φ1 |Ω1 ,1+α1 + S3 |ˆ
1/2+α2 +C2 h2
u|Ω2 ,1+α2 ) ((S2 + S4 )|φ2 |Ω2 ,1+α2 + S3 |ˆ (k)
with α1 /2 ≤ β ≤ α1 , S1 is a stability factor independent of h, S2 = kA1 ∂n U1 kΓ , S3 = kA1 ∂n φ1 kΓ , and S4 = kσ (k) kΓ .
Corollary 5.6. Assume the triangulations T1,h and T2,h match along the (k) (k) interface Γ. If U1 and U2 solve (3.3) and (3.4) respectively, then the errors
CONJUGATE HEAT TRANSFER (k)
e1 = u1 − U1
(k)
and e2 = u2 − U2
ke1 kΩ1 + ke2 kΩ2 ≤
satisfy
! 1 Ch1+α |φ1 |K,1+α1 K · 1 Ch1+α |φ1 |K,1+α1 K ! (k) 2 kf2 − L2 U2 kK Ch1+α |φ2 |K,1+α2 K · −1/2 (k) 2 1 Ch1+α |φ2 |K,1+α2 [A2 ∂n U2 ]k∂K K 2 khK −1/2 2 |φ2 |1+α2 + kh2 [A∂n U (k) ]kΓ · Ch1+α 2 (k−1) (k) + (kA1 ∂n φ1 kΓ ) · kU2 − U2 kΓ 1 kf1 kΩ1 · (kπ2 φ2 kΓ ) , + Chβ1 kf1 kΩ1 + CS1 h1+α 1 (k)
X
K∈T1,h
+
17
X
K∈T2,h
kf1 − L1 U1 kK −1/2 (k) 1 kh [A1 ∂n U1 ]k∂K K 2
with α1 /2 ≤ β ≤ α1 and S1 is a stability factor independent of h1 . It is clear that the term containing hβ1 decreases at a slower rate than the other terms. In fact, if we assume u1 ∈ H 2 (Ω1 ) and u2 ∈ H 2 (Ω2 ) then the other terms are O(h21 ) or O(h22 ), while the finite element flux is O(h1 ) at best. Suppose that instead we solve (3.8) rather than (3.4), i.e., we pass σ (k) instead of the finite element flux. This changes the fourth term in the error representation formula to I4 = (σ (k) , π2 φ2 )Γ − (σ (k) , π1 φ1 )Γ .
Theorem 5.7. Assume the triangulations T1,h and T2,h do not match along (k) (k) the interface Γ. If U1 and U2 solve (3.3) and (3.8) respectively, then the errors (k) (k) e1 = u1 − U1 and e2 = u2 − U2 satisfy ke1 kΩ1 + ke2 kΩ2 ≤
K∈T1,h
+
X
! 1 Ch1+α |φ1 |K,1+α1 K · 1 Ch1+α |φ1 |K,1+α1 K ! (k) 2 Ch1+α |φ2 |K,1+α2 kf2 − L2 U2 kK K · −1/2 (k) 2 1 Ch1+α |φ2 |K,1+α2 [A2 ∂n U2 ]k∂K K 2 khK −1/2 2 |φ2 |1+α2 + kh2 [A∂n U (k) ]kΓ · Ch1+α 2 (k−1) (k) + (kA1 ∂n φ1 kΓ ) · kU2 − U2 kΓ (k)
X
K∈T2,h
kf1 − L1 U1 kK −1/2 (k) 1 kh [A1 ∂n U1 ]k∂K K 2
1/2+α1
+C1 h1
1/2+α2 +C1 h2
u|Ω1 ,1+α1 ) ((S2 + S4 )|φ1 |Ω1 ,1+α1 + S3 |ˆ
u|Ω2 ,1+α2 ) , ((S2 + S4 )|φ2 |Ω2 ,1+α2 + S3 |ˆ
(k)
where S2 = kA1 ∂n U1 kΓ , S3 = kA1 ∂n φ1 kΓ , and S4 = kσ (k) kΓ . Corollary 5.8. Assume the triangulations T1,h and T2,h match along the (k) (k) interface Γ. If U1 and U2 solve (3.3) and (3.8) respectively, then the errors
18
D. ESTEP, S. TAVENER, AND T. WILDEY (k)
e1 = u1 − U1
(k)
and e2 = u2 − U2
ke1 kΩ1 + ke2 kΩ2 ≤
X
K∈T1,h
+
X
K∈T2,h
satisfy
! 1 Ch1+α |φ1 |K,1+α1 K · 1 Ch1+α |φ1 |K,1+α1 K ! (k) 2 Ch1+α |φ2 |K,1+α2 kf2 − L2 U2 kK K · −1/2 (k) 2 1 Ch1+α |φ2 |K,1+α2 [A2 ∂n U2 ]k∂K K 2 khK −1/2 2 |φ2 |1+α2 + kh2 [A∂n U (k) ]kΓ · Ch1+α 2 (k−1) (k) + (kA1 ∂n φ1 kΓ ) · kU2 − U2 kΓ . (k)
kf1 − L1 U1 kK −1/2 (k) 1 kh [A1 ∂n U1 ]k∂K K 2
Comparing Theorem 5.7 with Theorem 5.5 and Corollary 5.8 with Corollary 5.6, we see that the terms involving hβ have dropped out and the optimal order of convergence has been restored. 5.2. Numerical results. Example 5.1. Let Ω1 = [0, 1] × [0, 1] and the triangulations match along Γ = {(x, y)| x = interface problem −∇ ( · (∇u1 ) = f1 , u =u , 1 2 ∂ u = ∂n u 2 , n 1 −∇ · (∇u2 ) = f2 ,
Ω2 = [1, 2] × [0, 1] and assume that 1, 0 ≤ y ≤ 1}. Consider the elliptic x ∈ Ω1 , x ∈ Γ,
(5.2)
x ∈ Ω2 ,
where the data f1 (x, y) and f2 (x, y) are chosen so the true solutions are u1 = u2 = sin(πx/2) sin(2πy). The purpose of this experiment is to show how the asymptotic O(hβ1 ) convergence rate may be masked by apparently minor changes to the problem, and therefore may be difficult to observe in all situations. We change the boundary conditions along the top boundary, solving (5.2) with the boundary conditions u1 = 0, along x = 0 and y = 0, ∂ u = 2π sin(πx/2), along y = 1, n 1 u along x = 2 and y = 0, 2 = 0, ∂n u2 = 2π sin(πx/2), along y = 1. We set α = 1/2 and iterate until the convergence criteria is met. For comparison, we also solve (5.2) using a fully coupled method on the same mesh in COMSOL. We see, in Figure 5.1, that the fully coupled solution converges quadratically, while the iterative solution converges linearly. Next, we solve (5.2) with Dirichlet boundary conditions ( u1 = 0, x ∈ ∂Ω1 \Γ, u2 = 0, x ∈ ∂Ω2 \Γ, and again we also solve the fully coupled problem for comparison. In Figure 5.2, we observe that the operator decomposition solution initially converges quadratically, but eventually the O(hβ1 ) term dominates and the convergence is ultimately linear.
19
CONJUGATE HEAT TRANSFER
Fully Coupled Solution Operator Decomposition
Fully Coupled Solution Operator Decomposition
−3
−4
log( L2error)
log( L2 error)
−3
−5 −6
−4 −5 −6
−7 1.5
2 2.5 3 log(sqrt(degrees of freedom))
−7 1.5
3.5
2 2.5 3 log(sqrt(degrees of freedom))
3.5
Fig. 5.1. Comparison of L2 error in the fully coupled approximation and the operator decomposition approximation over Ω1 (left) and Ω2 (right) with mixed boundary conditions on ∂Ωi \Γ for i = 1, 2. −2
−2 Fully Coupled Solution Operator Decomposition
Fully Coupled Solution Operator Decomposition
−4 log(L2error)
log(L2 error)
−4 −6 −8 −10
−6 −8 −10
2
3
4
5
6
log(sqrt(degrees of freedom))
2
3 4 5 log(sqrt(degrees of freedom))
6
Fig. 5.2. Comparison of L2 error in the fully coupled approximation and the iterative approximation over Ω1 (left) and Ω2 (right) with Dirichlet boundary conditions on ∂Ωi \Γ for i = 1, 2.
Example 5.2. In the final example we assume the triangulations do not match along the interface and consider (5.2) with the mixed boundary conditions u1 = 0, A1 ∂n u1 A ∂ u 1 n 1 u = 0, 2 A2 ∂n u2 A2 ∂n u2
= 2π sin(πx/2), = −2π sin(πx/2), = 2π sin(πx/2), = −2π sin(πx/2),
along along along along along along
x = 0, y = 1, y = 0, x = 2, y = 1, y = 0. (k)
First, we solve the problem iteratively by passing the finite element flux, A1 ∂n U1 . Next, we use the boundary flux method to compute and pass σ (k) . In Figure 5.3, we compare the L2 errors over Ω1 ∪ Ω2 on a series of meshes. We see that the approximations using the finite element flux converge linearly, while the approximations using the boundary flux converge quadratically. 6. Conclusion. We have conducted an a posteriori error analysis of a finite element implementation of a number of operator decomposition strategies for a canonical
20
D. ESTEP, S. TAVENER, AND T. WILDEY
Using Finite Element Flux Using Boundary Flux
−5
2
log(L error)
−4
−6
−7 2.8
3
3.2 3.4 3.6 3.8 4 log(sqrt(degrees of freedom))
4.2
4.4
Fig. 5.3. Comparison of L2 error in the iterative approximation using the finite element flux and the boundary flux recovery method over Ω1 ∪ Ω2 with mixed boundary conditions on ∂Ωi \Γ for i = 1, 2.
conjugate heat transfer problem. By modifying the standard approach based on variational analysis, residuals and the generalized Green’s function, we derive accurate error estimates and use these estimates to guide adaptive mesh refinement. Modifications to the standard error analysis account for the transfer of error between components of the decomposed operator and interpolation error. Our analysis accounts for the differences between the adjoints of the full problem and the operator decomposition system. We show that the loss of order typically observed in the operator decomposition method is due to inaccuracy of transferred gradient information and we show how the boundary flux method can be used to efficiently regain the expected optimal order of convergence. REFERENCES [1] I. Babuska. The finite element method for elliptic equations with discontinuous coefficients. Computing, 5:207–213, 1970. [2] Ivo Babuska and Manil Suri. The p and hp versions of the finite element method, basic principles and properties. SIAM Review, 36:578–632, 1994. [3] Wolfgang Bangerth and Rolf Rannacher. Adaptive Finite Element Methods for Differential Equations. Birkhauser Verlag, 2003. [4] R. Becker and R. Rannacher. An optimal control approach to a posteriori error estimation in finite element methods. Acta Numerica, pages 1–102, 2001. [5] James H. Bramble and J. Thomas King. A finite element method for interface problems in domains with smooth boundaries and interfaces. Advances in Comp. Math., 67:1–19, 1998. [6] T. Cao, D.W. Kelly, and M. Ainsworth. Some useful techniques for pointwise and local error estimates of the quantities of interest in the finite element approximation. ANZIAM, 42:317–339, 2000. [7] G.F. Carey. Derivative calculation from finite element solutions. Comp. Meth. in Applied Mech. and Engr., 35:1–14, 1982. [8] V. Carey, D. Estep, and S. Tavener. A-posteriori error control of one-way coupled elliptic systems. In preparation, 2006. [9] Zhiming Chen and Jun Zou. Finite element methods and their convergence for elliptic and parabolic interface problems. Numer. Math, 79:175–202, 1998. [10] K. Eriksson, D. Estep, P. Hansbo, and C. Johnson. Introduction to adaptive methods for differential equations. Acta Numerica, pages 105–158, 1995. [11] K. Eriksson, D. Estep, P. Hansbo, and C. Johnson. Computational Differential Equations. Cambridge University Press, New York, 1996. [12] D. Estep. A posteriori error bounds and global error control for approximations of ordinary differential equations. SIAM Journal on Numerical Analysis, 32:1–48, 1995. [13] Donald Estep. A short course on duality, adjoint operators, green’s functions, and a-posteriori
CONJUGATE HEAT TRANSFER
21
error analysis. Sandia National Laboratories, Albuquerque, New Mexico, 2004. [14] Donald Estep, Michael Holst, and Mats Larson. Generalized Green’s functions and the effective domain of influence. SIAM Jour. Sci. Comp., 2004. [15] Donald Estep, Roy Williams, and Mats Larson. Estimating the error of numerical solutions of systems of reaction-diffusion equations. Memoirs of the American Mathematical Society, 146, 2000. [16] M.K. Seager G.F. Carey, S.S. Chow. Approximate boundary-flux calculations. Comp. Meth. in Applied Mech. and Engr., 50:107–120, 1985. [17] M. Giles, M. G. Larson, J. M. Levenstam, and E. Suli. Adaptive error control for finite element approximations of the lift and drag coefficients in viscous flow. Technical Report NA-97-06, Oxford University Computing Laboratory, 1997. [18] M.B. Giles. Stability analysis of numerical interface boundary conditions in fluid-structure thermal analysis. International Journal for Numerical Methods in Fluids, 25:421–436, 1997. [19] V. Ginting, D. Estep, J. Shadid, and S. Tavener. A-posteriori analysis of operator splitting for ordinary differential equations. In preparation, 2006. [20] V. Heuveline and R. Rannacher. Duality-based adaptivity in the hp-finite element method. J. Numer. Math, 0:1–18, 2003. [21] Cornelius Lanczos. Linear Differential Operators. Dover Publications, 1997. [22] Ahmed Naga and Zhimin Zhang. A posteriori error estimates based on the polynomial preserving recovery. SIAM J. Numer. Anal., 42:1780–1800, 2004. [23] J.R. Rice, P. Tsompanopoulus, and E. Vavalis. Interface relaxation methods for elliptic differential equations. Applied Numerical Mathematics, 32:219–245, 2000. [24] Barry Smith, Petter Bjorstad, and William Gropp. Domain Decomposition. Cambridge University Press, 1994. [25] M.F. Wheeler. A Galerkin procedure for estimating the flux for two-point boundary-value problems using continuous picewise-polynomial spaces. Numer. Math, 2:99–109, 1974. [26] T. Wildey. A posteriori analysis of operator decomposition in interface problems. Ph.D. Thesis, Colorado State University, 2007. [27] T. Wildey, D. Estep, and S. Tavener. A posteriori error estimation of approximate boundary fluxes. In preparation, 2006. [28] Daoqi Yang. A parallel nonoverlapping schwarz domain decomposition method for elliptic interface problems. IMA Journal on Numerical Analysis, 16:75–91, 1996. [29] Ivan Yotov. A multilevel newton-krylov interface solver for multiphysics coupling of flow in porous media. Numer. Linear Algebra Appl., 8:551–570, 2001. [30] O.C. Zienkiewicz and J.Z. Zhu. The superconvergent patch recovery and a posteriori error estimates. Intert. J. Numer. Methods Engrg., 33:1331–1364 and 1365–1382, 1992.
