Study on Error Estimates and Adaptivity

The goal-oriented error estimator developed in [197, 199] is used as a refinement criterion. The error in the value of a linear target functional is measured in terms of weighted residuals that depend on the solutions to the primal and dual problems. The Galerkin orthogonality error is taken into account and turns out to be dominant whenever flux or slope limiters are activated to enforce monotonicity constraints. The localization of global errors is performed using a natural decomposition of the involved weights into nodal contributions. The developed simulation tools are applied to a linear convection problem in two space dimensions. The goal-oriented approach to error estimation [14, 27, 185, 295, 309] is applicable not only to elliptic PDEs but also to hyperbolic conservation laws [141, 142, 310]. In most cases, the error in the quantity of interest is estimated using the duality argument, Galerkin orthogonality, and a direct decomposition of the weighted residual into element contributions. The most prominent representative of such error estimators is the Dual Weighted Residual (DWR) method of Becker and Rannacher [27, 28].

The recent paper by Meidner et al. [248] is a rare example of a DWR estimate that does not require Galerkin orthogonality or information about the cause of its possible violation. Kuzmin and Korotov [197] applied the DWR method to steady convection diffusion equations and obtained a simple estimate of local Galerkin orthogonality errors due to flux limiting or other ‘variational crimes.’ In contrast to the usual approach, the weighted residuals are decomposed into nodal (rather than element) contributions. In regions of insufficient mesh resolution, the computable Galerkin orthogonality error comes into prominence. The mesh adaptation strategy to be presented below takes advantage of this fact. Steady convective transport of a conserved scalar quantity u boundary vu) = s Here v is a stationary velocity field and s is a volumetric source/sink. Due to hyperbolicity, a Dirichlet boundary condition is imposed at the inlet u = uD in = {x v ·n < 0}, (5.2) where n is the unit outward normal and uD is the prescribed boundary data. The weak form of the above boundary value problem can be written as a(w,u) = b(w), 8w. (5.3) For brevity, we refrain from an explicit definition of functional spaces. The bilinear form a(·, ·) and the linear functional b(·) are defined by The inflow boundary conditions are imposed weakly via the surface integrals.

Available online at www.ignited.in Page 2

The differentiation of vu in (5.4) can be avoided using integration by parts this representation implies that a discontinuous weak solution u is admissible. In linear hyperbolic problems of the form (5.1), singularities travel along the streamlines of v. They may be caused by a jump in the value of s or uD.

Let uh be a continuous function that may represent an approximate solution to (5.1)– (5.2) or a finite element interpolant of discrete nodal values. The numerical error e = u−uh can be measured using the residual of (5.3) w,uh) = b(w)−a(w,uh). (5.7) w,uh) depends not only on the quality of uh but also on the choice of w. In goal-oriented estimates, this weight carries information about the 5.3 Global Error Estimates 199 quantities of interest. The objectives of a numerical study are commonly defined in terms of a linear output functional, such as [310] The piecewise-constant function g picks out a subdomain, for example, an interior or boundary layer, where a particularly accurate approximation to u is desired. The selector h picks out a portion of the out= {x v · n > 0}, where the convective flux is to be controlled. In order to estimate the error j(e) in the numerical value of the output functional, consider the dual or adjoint problem [27, 28] associated with (5.3) a(z, e) = j(e), 8e. (5.9) The surface integral in (5.8) implies the weakly imposed Dirichlet boundary condition z = h j(e) and residual (5.7) are related by j(u−uh) = a(z,u−uh) z,uh). (5.10) An arbitrary numerical approximation zh to the exact solution z of the dual problem (5.9) can be used to decompose the so-defined error as follows j(u−uh) z−zh,uhzh,uh). (5.11) If Galerkin orthogonality holds for the numerical approximation uhzh,uh)= 0. Thus, the computable zh,uh) is omitted in most goal-oriented error estimates for finite element discretizations. However, the orthogonality condition is frequently violated due to numerical integration, round-off errors, slack tolerances for iterative solvers, and flux limiting. Since the exact dual solution z is usually unknown, the derivation of a computable error estimate involves another approximation ˆz _ z such that j(u−uh) z−zh,uhzh,uh). (5.12) The magnitudes of the two residuals can be estimated z−zh,uhzh,uh(5.13) assembled from contributions of individual nodes or elements, as explained in the next section. The reference solution ˆz is commonly obtained from zh zh,uh) = 0, then the estimate j(u−uh) _ 0 that follows from (5.12) with ˆz = zh is worthless, hence the need to compute ˆz on another mesh or interpolate it using higher-order polynomials [197, 295]. On the other hand, the setting ˆz =zh is not only acceptable but also optimal for nonlinear flux-limited discretizations such that j(u−uhzh,uh) 6= 0. In situations when the tz−zh,uh) is non200 5 Error Estimates and Adaptivity negligible, extra work needs to be invested into the recovery of a superconvergent approximation ˆz 6= zh.

verify the accuracy of the approximate solution uh but the estimated errors in the quantity of interest must be localized to find the regions where a given mesh is too coarse or too fine. A straightforward decomposition of weighted residuals into element contributions results in an oscillatory distribution and a strong overestimation of local zh,uh) to a k can be large in magnitude even if Galerkin orthogonality is satisfied globally (positive and negative contributions cancel out). Following Schmich and Vexler [295]

In this section, the presented high-resolution finite element scheme, goal-oriented error estimator, and hierarchical mesh adaptation algorithm are applied to a test problem from [156]. Consider equation (5.1) with s _ 0 and v(x, y) = (y,−x This incompressible velocity field corresponds to steady rotation about (0,0). The exact solution and inflow boundary conditions are given by [156] u(x, y) =_ 1, if 0.35 _ p x2+y2 _ 0.65,

Available online at www.ignited.in Page 3

0, otherwise. The so-defined discontinuous inflow profile (−1 _ x < 0, y = 0) undergoes circular convection and propagates along the streamlines of v(x, y) all the way to the outlet (0 < x _ 1, y = 0), while its shape remains the same. Let j(u) be defined by (5.8) with g (−0.1,0.1)×(0,1) and g = 0 elsewhere. The function h is defined as the trace of g j(u) is 6.04497e−02. The solution shown in Fig. 5.1 (a) was computed by the FEM-LED scheme described in Chapter 4 on a uniform mesh of bilinear elements with spacing h = 1/80. Owing to algebraic flux correction, the resolution of the discontinuous front is remarkably sharp, and no undershoots or overshoots are observed. However, it is obvious that there is actually no need for such a high resolution beyond x > 0.1 if it is enough to have an on the solution in this subdomain. This is illustrated by Fig. 5.1 (b) which shows the solution to the dual problem computed by the FEM-LED scheme on the same mesh.

limiting. Remarkably, the resulting global estimates are in a good agreement with the exact error which is illustrated in Table 5.1 for different grid spacings. The sharpness of the obtained error estimates is measured using the absolute and relative effectivity indices [197] We remark that the value of Ieff is unstable and misleading when the denominator is very small or zero, and the evaluation of integrals is subject to rounding errors. The relative effectivity index Ieff is free of this drawback and exhibits monotone convergence as the mesh is refined 202 5 Error Estimates and Adaptivity

The adaptive hybrid mesh presented in Fig. 5.2 is refined along the discontinuity lines of u but only until they cross line x = 0.1 would not improve the accuracy of the solution uh h = 1/320, which corresponds to more than 200,000 cells in the case of global mesh refinement. Since the dual weight zh contains built-in information regarding the transport of errors and goals of simulation, such error estimators furnish a better refinement criterion than, for example, error indicators based on gradient recovery [362]. In the latter case, unnecessary mesh refinement would take place along the discontinuities located downs

1. M. Aftosmis and N. Kroll, A quadrilateral based second-order TVD method for unstructured adaptive methods. AIAA Paper, 91-0124, 1991. 2. M. Ainsworth and J.T. Oden, A Posteriori Error Estimation in Finite Element Analysis. John Wiley & Sons, New York, 2000. 3. M. Ainsworth, J.Z. Zhu, A.W. Craig, O.C. Zienkiewicz, Analysis of the Zienkiewicz-Zhu aposteriori error estimator in the finite element method. Int. J. Numer. Methods Engrg. 28:9 (1989) 2161–2174. 4. J.D. Anderson, Jr., Computational Fluid Dynamics. The Basics with Applications. McGraw- Hill, 1995. 5. J.D. Anderson, Jr., Modern Compressible Flow: With Historical Perspective, McGraw-Hill, 1990. 6. F. Angrand, A. Derieux, L. Loth, G. Vijayasundaram, Simulation of Euler transonic flows by means of explicit finite element type schemes. INRIA Research Report 250, 1983. 7. F. Angrand and A. Derieux, Some explicit triangular finite element schemes for the Euler equations. Int. J. Numer. Methods Fluids 4 (1984) 749–764. 8. P. Arminjon and A. Dervieux, Construction of TVD-like artificial viscosities on 2- dimensional arbitrary FEM grids. INRIA Research Report 1111, 1989. 9. M. Arora and P.L. Roe, A well-behaved TVD limiter for high-resolution calculations of unsteady flow. J. Comput. Phys. 132 (1997) 3–11. 10. Athena test suite, http://www.astro.virginia.edu/VITA/ATHENA/dmr.html. 11. K. Baba and M. Tabata, On a conservative upwind finite element scheme for convective diffusion equations. RAIRO Numerical Analysis 15 (1981) 3–25. 12. I. Babuˇska and W.C. Rheinboldt, Error estimates for adaptive finite element computations. SIAM J. Numer. Anal. 15:4 (1978) 736–354.

Available online at www.ignited.in Page 4

13. I. Babuˇska and W.C. Rheinboldt, A posteriori error estimates for the finite element method. Int. J. Numer. Methods Engrg. 12:10 (1978) 1597–1615. 14. W. Bangerth and R. Rannacher, Adaptive finite element methods for differential equations. Lectures in Mathematics, ETH Z¨urich, Birkh¨auser, 2003. 15. R.E. Bank, A.H. Sherman, A. Weiser, Some refinement algorithms and data structures for regular local mesh refinement. In: R. Stepleman (ed.), Scientific Computing, Applications of Mathematics and Computing to the Physical Sciences, IMACS Transactions on Scientific Computation, Vol. I. North-Holland, Amsterdam, 1983, 3–17.

17. T.J. Barth, Numerical aspects of computing viscous high Reynolds number flows on unstructured meshes. AIAA Paper, 91-0721, 1991. 18. T.J. Barth, Aspects of unstructured grids and finite volume solvers for the Euler and Navier- Stokes equations. In: Lecture Series 1994-05, von Karman Institute for Fluid Dynamics, Brussels, 1994.