T-IFISS: a toolbox for adaptive FEM computation
TT-IFISS: a toolbox for adaptive FEM computation ∗ Alex Bespalov † Leonardo Rocchi † David Silvester ‡ Abstract
T-IFISS is a finite element software package for studying finite element solutionalgorithms for deterministic and parametric elliptic partial differential equations.The emphasis is on self-adaptive algorithms with rigorous error control using a va-riety of a posteriori error estimation techniques. The open-source MATLAB frame-work provides a computational laboratory for experimentation and exploration, en-abling users to quickly develop new discretizations and test alternative algorithms.The package is also valuable as a teaching tool for students who want to learn aboutstate-of-the-art finite element methodology.
Key words : finite elements, stochastic Galerkin methods, a posteriori error estimation,adaptive methods, goal-oriented adaptivity, PDEs with random data, parametric PDEs,mathematical software
AMS Subject Classification : 97N80, 65N30, 65N15, 65N50, 35R60, 65C20, 65N22
Progress in computational mathematics is frequently motivated and supported by theresults of numerical experiments. The well-established IFISS software package [26, 24]is associated with the monograph [25] and is structured as a stand-alone package forstudying discretization algorithms for partial differential equations (PDEs) arising in in-compressible fluid dynamics. IFISS is also an established starting point for developingcode for specialized research applications. The package is currently used in universitiesaround the world to enhance the teaching of advanced courses in mathematics, compu-tational science and engineering. Investigative numerical experiments enable students todevelop deduction and interpretation skills and are especially useful in helping studentsto remember critical ideas in the long term.The
T-IFISS (Triangular IFISS) toolbox extends the IFISS philosophy and design fea-tures to finite element approximations on triangular grids for deterministic and stochas-tic/parametric elliptic partial differential equations. The emphasis of T-IFISS is on self-adaptive algorithms with rigorous error control using a variety of a posteriori error es-timation techniques. In the same way as for its precursor, the development of T-IFISS ∗ This work was supported by the EPSRC under grants EP/P013317/1 and EP/P013791/1, and waspartially supported by The Alan Turing Institute under the EPSRC grant EP/N510129/1. † School of Mathematics, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK( [email protected] , [email protected] ). ‡ Department of Mathematics, University of Manchester, Oxford Road, Manchester M13 9PL, UK( [email protected] ). IFISS is an acronym for “Incompressible Flow & Iterative Solver Software”. See the swMATH resource page http://swmath.org/software/4398 . a r X i v : . [ m a t h . NA ] M a r as been motivated by a desire to create a computational laboratory for experimentationand a tool for reproducible research in computational science; indeed, the toolbox enablesusers not only to quickly explore new discretization strategies and test alternative algo-rithms, but to replicate, validate and independently verify computational results. Thus,in this article, instead of giving a comprehensive technical description of the software,we aim to highlight and document those features of T-IFISS that are not available in itsprecursor package IFISS; we will illustrate these features with a series of case studies thatdemonstrate the efficiency of the software and its utility as a research and teaching tool.All the test problems that are built into the current version of T-IFISS are associ-ated with steady-state diffusion equations with anisotropic (or uncertain) conductivitycoefficients together with mixed (Dirichlet and Neumann) boundary conditions. ThesePDE problems may be solved on general polygonal domains, including slit domains anddomains with holes. It is worth pointing out that many of these (deterministic) testproblems could also be solved using general purpose finite element software packages likedeal.II [4], DUNE [17], FEniCS [39], and FreeFEM [34]. The attraction of the T-IFISSenvironment is the ease with which one can test the alternative error estimation, marking,and refinement strategies, develop new strategies, and extend functionality to new prob-lem classes as well as new approximations and solution algorithms. The modular structureof the code allows one to examine interactions between different components in order toachieve favourable convergence properties or minimize the associated computational cost.Finally, the use of the high-level MATLAB syntax ensures readability, accessibility, andquick visualization of the solution, error estimators, and finite element meshes. Thesefeatures are invaluable in teaching but are not so explicit in other adaptive finite elementmethod (FEM) packages like ALBERTA [51], PLTMG [5], and p1afem [29]. A uniquefeature of T-IFISS is that it can be used to solve parameter-dependent elliptic PDE prob-lems stemming from uncertainty quantification models. This facility, called stochastic T-IFISS [14], develops the idea of adaptive stochastic Galerkin FEM and provides aneffective alternative to traditional sampling methods commonly used for such problems.The only software packages that we know of that have a similar capability are ALEA [23]and SGLib [56].The following is a brief overview of the functionality available in T-IFISS (includingStochastic T-IFISS) with links to its directory structure (we refer to [50, Appendix B] formore details including the description of the associated data structures): • computing Galerkin solutions with P and P approximations over the specifiedtriangular mesh (directory /diffusion ); • computing stochastic Galerkin FEM solutions with spatial P and P approxima-tions over the specified triangular mesh and for the specified polynomial approxi-mation on the parameter domain; the functionality here includes effective iterativesolution of very large linear systems stemming from stochastic Galerkin FEM (di-rectory /stoch diffusion ); • estimating the energy error in the computed Galerkin solution and the error in aquantity of interest derived from this solution (directories /diffusion error and /goafem , respectively); T-IFISS version 1.2 was released in February 2019 and runs under MATLAB or Octave. It canbe downloaded from and is compatible with Win-dows, Linux and MacOS computers. estimating the energy error in the computed stochastic Galerkin solution and the er-ror in a quantity of interest derived from this solution (directories /stoch diffusion and /stoch diffusion/stoch goafem , respectively); • adaptive refinement of Galerkin approximations, including local adaptive mesh re-finement and, in the case of stochastic/parametric problems, adaptive enrichmentof multivariable polynomial approximations (directories /diffusion adaptive and /stoch diffusion/stoch diffusion adapt ); • visualization of Galerkin solutions, goal functionals, error estimators, finite elementmeshes; plotting convergence history (the error estimates against the number ofdegrees of freedom) (directory /graphs );The rest of the article is organized as follows. In the next section, we recall the mainingredients of adaptive FEM and describe their implementation in T-IFISS; we illustratethis with two case studies (solving the diffusion equation with strongly anisotropic coeffi-cient and computing a harmonic function in the L-shaped domain). Section 3 is focusedon implementing the goal-oriented error estimation and adaptivity; the case study heredemonstrates the capability of the software to approximate the value of a singular solutionto the diffusion problem at a fixed point away from the singularity. In section 4, we discussthe implementation of stochastic Galerkin FEM including the associated error estimationand adaptivity. The efficiency of our adaptive algorithm is demonstrated with a repre-sentative case study (solving the diffusion equation with parametric coefficient over theL-shaped domain). Some extensions and future developments of the toolbox are brieflyoutlined in Section 5. Adaptive finite element approximations to solutions of partial differential equations aretypically computed by iterating the following loop of four component modules:
SOLVE = ⇒ ESTIMATE = ⇒ MARK = ⇒ REFINE . (2.1)In this section we first describe the implementation of the four components in (2.1)within T-IFISS. This is followed by two case studies that highlight some of the featuresof the toolbox and demonstrate its utility and its efficiency. Module
SOLVE . For a given PDE problem, the Galerkin solution defined on a specific(structured or unstructured ) triangular mesh is computed by solving the linear systemassociated with the Galerkin projection of the variational formulation of the problem ontothe corresponding finite element space. Two types of ( C ) finite element approximationsare implemented: piecewise linear ( P ) and piecewise quadratic ( P ). For either of theseapproximations, fast computation of the entries in the stiffness matrix and the load vector In the MATLAB release of T-IFISS, unstructured meshes are generated using the DistMesh pack-age [46].
3s achieved by vectorizing the calculations over all the elements. When solving a deter-ministic problem the resulting linear equation system is solved using the highly optimizedsparse direct solver (UMFPACK) that is built into MATLAB and Octave. Module
ESTIMATE . The purpose of this module is two-fold. First, it computes localerror indicators that provide information about the distribution of estimated local errorsin the computed Galerkin solution; the error indicators may be associated with elementsor edges of the underlying triangulation. Second, it computes an estimate of the (total)energy error in the Galerkin solution. This estimate is used to decide whether the stoppingtolerance is met. T-IFISS offers a choice of the following three error estimation strategies(EES) for linear ( P ) approximation.(EES1) is a local hierarchical error estimator computed via a standard element residualtechnique (see [1, Section 3.3]) using either piecewise linear or piecewise quadratic bubblefunctions over subelements obtained by uniform refinements . The local error indicatorsin this case are computed elementwise by solving 3 × (cid:96) -norm of thevector of local error indicators.(EES2) is a global hierarchical estimator (see [6], [1, Section 5]) using piecewise linearbubble functions corresponding to the uniform refinement of the original triangulation.Note that the implementation of this strategy requires solving a sparse linear systemassociated with a global residual problem. The localizations of the estimator (to eitherthe elements (default option) or the edges of triangulation) gives two types of local errorindicators in this EES.(EES3) is a two-level error estimate employing piecewise linear bubble functions asso-ciated with edge midpoints (of the original triangulation); see [44, 43]. In this case, it isnatural to choose local error indicators associated with edges (this is the default choice in(EES3)). However, the choice of elementwise error indicators is also offered as an option;these are computed for each interior element from three corresponding edge indicators (orfrom two indicators for the elements with an edge on the Dirichlet part of the boundary).There is currently no flexibility with choosing the EES when using quadratic ( P )approximation: the local hierarchical error estimation strategy (EES1) is employed withpiecewise quartic bubble functions. More specifically, the local error indicators are com-puted elementwise by solving 9 × (cid:96) -norm of the vector oflocal error indicators.Module MARK . In this module the elements (or edges) with largest error indicatorsare selected (i.e., marked) for refinement. Two popular marking strategies are currentlyimplemented: the maximum strategy and the
D¨orfler strategy (also referred to as theequilibration or bulk chasing strategy).Let { β ( s ); s ∈ S} denote the set of error indicators associated with the elementsof the set S (e.g., S can be the set of edges or elements of the triangulation). In themaximum marking strategy, that dates back to [3], the element s ∈ S is marked if the We would almost certainly use an iterative solver preconditioned with an algebraic multigrid V-cycleif we were trying to solve the same PDE problem in three spatial dimensions. The default uniform refinement in T-IFISS is by three bisections (see Figure 1(d)). However, thereis an option to switch to the so-called red uniform refinement (i.e., the one obtained by connectingthe edge midpoints of each triangle); this can be done by setting subdivPar = 1 within the function adiff adaptive main.m . More precisely, the interior edges and the edges associated with those parts of the boundary wherethe Neumann and non-homogeneous Dirichlet boundary conditions are prescribed. a) (b) (c) (d) Figure 1. NVB bisections: (a) one, (b)-(c) two, and (d) three bisections of the edges of the triangularelement. Double lines indicate reference edges, black dots indicate the edges to be bisected, and whitedots indicate the newest vertices. associated error indicator β ( s ) is larger than a fixed proportion of the maximum amongall error indicators. Specifically, for a given marking (or, threshold) parameter θ ∈ [0 , M ⊆ S of marked elements such that β ( s ) ≥ θ max s ∈S β ( s ) ∀ s ∈ M . (2.2)Note that in this strategy, smaller values of θ lead to larger subsets M .In the D¨orfler marking strategy, that was originally introduced in [20], sufficientlymany elements are marked such that the combined contribution of the corresponding errorindicators is larger than a fixed proportion of the total error estimate. More precisely,given a marking parameter θ ∈ (0 , M ⊆ S of minimalcardinality such that { β ( s ); s ∈ M} is the set of M largest error indicators and (cid:88) s ∈M β ( s ) ≥ θ (cid:88) s ∈S β ( s ) . (2.3)Here, smaller values of θ lead to smaller subsets M .Module REFINE . Given the set of marked elements (or, marked edges) that is obtainedby employing one of the above marking strategies, local adaptive mesh refinement isperformed in T-IFISS by implementing the longest edge bisection (LEB) strategy—avariant of the newest vertex bisection (NVB) method (we refer to [52, 49, 7, 38, 54, 45]for theoretical and implementational aspects of NVB refinements, as well as to [41] foran overview and comparison of NVB with other mesh-refinement techniques). In thismethod, a reference edge is designated for each triangle T (for the coarsest mesh, thisis always the longest edge of T ), and T is bisected by halving the reference edge; seeFigure 1(a). This introduces two new elements, the sons of T , for which reference edgesare selected . A recursive application of this procedure leads to a conforming mesh, whereone, two, or three bisections of the triangle T may be performed; see Figure 1(a)–(d). Therefinement by three bisections (see Figure 1(d)) is called the bisec3 refinement. Refiningall elements of the given mesh by three bisections results in a (conforming) uniform bisec3 refinement of this mesh.It is important to emphasize that NVB iteratively refines individual elements by bi-secting some (or all) of their edges. Therefore, either the set of marked elements or the In the NVB method, the reference edges of the sons are the edges opposite to the new vertex, whereasin the LEB method, the reference edge is always the longest edge of the element. Note, however, thatfor structured triangulations of square, L-shaped, and crack domains, both methods result in identicalrefinement patterns.
MARK , the mesh-refinement routine identifies the sub-sets of elements where one, two, or three bisections should be performed (see Figure 1)and then the elements in the three separate subsets are refined simultaneously.
Let us demonstrate the performance of adaptive finite element routines in T-IFISS withtwo test examples. When doing this, we will illustrate some of the ingredients of adaptiveFEM described in the previous subsection.
Example 1.
Let D = ( − , be the square domain. We consider the diffusionequation with a strongly anisotropic coefficient, together with a constant source functionand a homogeneous Dirichlet boundary condition: −∇ · ( A ∇ u ( x )) = 1 , x = ( x , x ) ∈ D,u ( x ) = 0 , x ∈ ∂D, (2.4)where A = [ ]. We solve this problem using P approximation. The solution isdepicted in Figure 2(b) and exhibits sharp gradients within the boundary layers along theedges x = ± tol = 1e-3 andemployed the element-based D¨orfler marking strategy (2.3) with θ = 0 .
5. The results ofcomputations are presented in Figures 3, 4 and in Table 1.Figure 3 shows convergence plots for the energy error estimates computed via eachof the three strategies. In Figure 4, we plot the locally refined meshes generated by − − T = 128 (a) − − · − (b)Figure 2. Example 1: (a) initial coarse mesh; (b) Galerkin solution to problem (2.4). − − − O ( N − . )degrees of freedom e rr o r e s t i m a t e (EES1) (a) O ( N − . )degrees of freedom (EES2) (b) O ( N − . )degrees of freedom (EES3) (c)Figure 3. Example 1: error estimates at each iteration of the adaptive algorithm employing the errorestimation strategies (EES1)–(EES3) and using the element-based D¨orfler marking with θ = 0 .
5. Here, N denotes the number of degrees of freedom. the adaptive algorithm in each of the three cases (for some intermediate tolerance). InTable 1, we have collected the data on iteration counts and mesh refinements for eachrun of the algorithm. It is evident from these results that using (EES1) leads to unstablereductions in the estimated errors and, as a consequence, to a large number of iterationsand an over-refined final mesh. This is because for problem (2.4) with constant coefficients,the elementwise interior residuals for P approximations have zero contributions to theassociated error estimator. In contrast, the adaptive algorithms employing (EES2) and(EES3) lead to essentially monotonic decay of the error estimates; the number of iterationsneeded to reach the stopping tolerance is nearly the same and similar mesh refinementpatterns are generated in both cases. We note that for all three strategies, the errorestimates decay with an optimal rate of O ( N − . ), where N is the number of degrees offreedom (d.o.f.).In our second experiment, we ran the adaptive algorithm driven by two-level errorestimates (i.e., using (EES3)) and employed the edge-based D¨orfler marking (2.3) with θ = T = 3200 (a) T = 1464 (b) T = 1462 (c)Figure 4. Example 1: locally refined meshes produced by the adaptive algorithm employing the errorestimation strategies (EES1)–(EES3) and using the element-based D¨orfler marking with θ = 0 .
5. Theheader T (cid:96) refers to the number of elements in the mesh at step (cid:96) of the adaptive process. daptive FEM with element-based D¨orfler marking ( θ = 0 . tol=1e-3 (EES1) (EES2) (EES3) L
77 25 24 η L . . . T L n L L denotesthe total number of iterations, η L , T L , and n L refer to the final error estimate, the number of elementsin the final mesh, and the number of degrees of freedom at the final iteration, respectively. . P approximation with a tightertolerance of tol = 2e-5 ). In Figure 5(a), we plot the energy error estimates at eachiteration. Comparing this plot with the one in Figure 3(c), we can see improvements interms of the monotonicity of the error decay and in terms of the number of iterations. Thecomputed effectivity indices for each iteration of the algorithm are plotted in Figure 5(b). − − O ( N − . )degrees of freedom e rr o r e s t i m a t e (EES3) (a) . . . O ( N − . )degrees of freedom e ff ec t i v i t y i nd e x (EES3) (b)Figure 5. Example 1: (a) error estimates at each iteration of the adaptive algorithm employing (EES3)and the edge-based D¨orfler marking with θ = 0 .
5, (b) the associated effectivity indices.
In our third experiment, we repeated the previous adaptive procedure with two smallerstopping tolerances. The overall computational times together with the contributions foreach of the components in (2.1) are reported in Table 2. Example 2.
This experiment addresses the question posed by Nick Trefethen and AbiGopal to the readers of the NA Digest in November 2018 (see [55, 33]). The communitywas challenged to compute (to a high accuracy) the point-value close to singularity fora harmonic function in the L-shaped domain. More precisely, let D = ( − , \ ( − , All timings reported in this paper were recorded using MATLAB 9.2 (R2017a) on a laptop with IntelCore i7 2.9GHz CPU and 16GB of RAM. daptive FEM with edge-based D¨orfler marking ( θ = 0 . tol=1e-3 tol=5e-4 tol=1e-4 L
16 19 27 η L . . . n L t ( SOLVE ) 0 .
167 0 .
811 49 . t ( ESTIMATE ) 0 .
324 2 .
013 107 . t ( MARK ) 0 .
004 0 .
016 0 . t ( REFINE ) 0 .
101 0 .
366 15 . t (overall) 3 .
090 7 .
937 439 . L , η L , and n L are as in Table 1. All times t are inseconds and the timings for individual modules are recorded at the final adaptive step. and consider the following problem: −∇ u ( x ) = 0 , x = ( x , x ) ∈ D,u ( x ) = (1 − x ) , x ∈ ∂D. (2.5)The goal is to compute u (0 . , .
01) to at least 8-digit accuracy (the exact value is1 . . . . ; see [33]).In this example, we set the stopping tolerance tol = 4e-5 and ran the adaptiveFEM algorithm with P approximation together with element-based D¨orfler marking with θ = 0 .
5. The prescribed tolerance was satisfied after 38 iterations (final number of d.o.f.was 253,231, run time 59.6 sec), giving the value u (0 . , . ≈ . − − (a) − − − − O ( N − )degrees of freedom e rr o r e s t i m a t e (b) − − T = 4862 (c)Figure 6. Example 2: (a) Galerkin solution to problem (2.5); (b) error estimates at each iteration of theadaptive algorithm with P approximations and tol = 4e-5 ; (c) adaptively refined mesh generated bythe algorithm for an intermediate tolerance of . Goal-oriented adaptivity
The error estimation strategies described in the previous section provide a mechanism forcontrolling approximation errors in the (global) energy norm. However, in many practicalapplications, simulations often target a specific quantity of interest (typically, a localfeature of the solution) called the goal functional . In such cases, the energy norm ofthe error is likely to be of limited interest. The implementation of goal-oriented errorestimation and adaptivity in T-IFISS is the focus of this section. We also discuss arepresentative case study to show the efficiency of the adopted approach.
We start by describing a general idea of the goal-oriented error estimation strategy im-plemented in T-IFISS. Let V be a Hilbert space and denote by V (cid:48) its dual space. Let B : V × V → R be a continuous, elliptic, and symmetric bilinear form with the associatedenergy norm |(cid:107) · (cid:107)| , i.e., |(cid:107) v (cid:107)| := B ( v, v ) for all v ∈ V . Given two continuous linearfunctionals F, G ∈ V (cid:48) , our aim is to approximate G ( u ), where u ∈ V is the unique solutionof the primal problem : B ( u, v ) = F ( v ) for all v ∈ V. (3.1)To this end, the standard approach (see, e.g., [27, 48, 9, 32]) considers z ∈ V as a uniquesolution to the dual problem : B ( v, z ) = G ( v ) for all v ∈ V . (3.2)Let V h be a finite dimensional subspace of V . Let u h ∈ V h (resp., z h ∈ V h ) be a uniqueGalerkin approximation of the solution to the primal (resp., dual) problem, i.e., B ( u h , v h ) = F ( v h ) (resp., B ( v h , z h ) = G ( v h )) for all v h ∈ V h . Then, it follows that | G ( u ) − G ( u h ) | = | B ( u − u h , z ) | = | B ( u − u h , z − z h ) | ≤ |(cid:107) u − u h (cid:107)| |(cid:107) z − z h (cid:107)| , (3.3)where the second equality holds due to Galerkin orthogonality.Assume that µ = µ ( u h ) and ζ = ζ ( z h ) are reliable estimates for the energy errors |(cid:107) u − u h (cid:107)| and |(cid:107) z − z h (cid:107)| , respectively, i.e., |(cid:107) u − u h (cid:107)| (cid:46) µ and |(cid:107) z − z h (cid:107)| (cid:46) ζ (here, a (cid:46) b means the existence of a generic positive constant C such that a ≤ Cb ).Hence, inequality (3.3) implies that the product µ ζ is a reliable error estimate for theapproximation error in the goal functional: | G ( u ) − G ( u h ) | (cid:46) µ ζ. (3.4) Having computed two Galerkin solutions u h , z h , the corresponding energy error estimates µ ( u h ), ζ ( z h ), and the reliable estimate µ ( u h ) ζ ( z h ) of the error in the goal functional(see (3.4)), a goal-oriented adaptive FEM (GOAFEM) algorithm proceeds by executingthe MARK and
REFINE modules of the standard adaptive loop (2.1).10hile the module
REFINE simply performs local mesh refinement as explained inthe previous section, the marking procedure in the GOAFEM algorithm requires specialcare. Specifically, since the error in the goal functional is controlled by the product ofthe energy error estimates for two Galerkin approximations (the primal and dual ones),the edge-marking for bisection (or, the element-marking for refinement) must take intoaccount the local error indicators associated with both approximations. Thus, given theset of local error indicators associated with the primal (resp., dual) Galerkin solution u h (resp., z h ), let M u (resp., M z ) denote the set of element edges that would be marked forbisection in order to enhance this Galerkin solution (to that end, one can use, e.g., theD¨orfler marking strategy, see (2.3)). There exist several strategies for combining the twosets M u and M z into a single marking set M that is used for mesh refinement in thegoal-oriented adaptive algorithm. Four different strategies are implemented in T-IFISS:(GO–MARK1) following [35], the marking set M is simply the union of M u and M z ;(GO–MARK2) following [42], the marking set M is defined as the set of minimalcardinality between M u and M z ;(GO–MARK3) following [8], the set M is obtained by performing D¨orfler marking onthe set of combined error indicators β ( E ) associated with edges of the triangulation, where β ( E ) := (cid:0) µ E ζ + ζ E µ (cid:1) / and µ E (resp., ζ E ) is the local contribution to µ (resp., ζ ) associated with the edge E ;(GO–MARK4) this marking strategy is a modification of (GO–MARK2); follow-ing [28], we compare the cardinality of M u and that of M z to define M (cid:63) := M u and M (cid:63) := M z if M u ≤ M z , M (cid:63) := M z and M (cid:63) := M u otherwise;the marking set M is then defined as the union of M (cid:63) and those M (cid:63) edges of M (cid:63) thathave the largest error indicators.Comparing these four strategies, it is proved in [28, Theorem 13] that the GOAFEMalgorithm employing marking strategies (GO–MARK2)–(GO–MARK4) generates approx-imations that converge with optimal algebraic rates, whereas only suboptimal convergencerates have been proved for marking strategy (GO–MARK1); cf. [28, Remark 4] and [35,Section 4]. The numerical results in [28] suggest that (GO–MARK4) is more effectivethan the original strategy (GO–MARK2) in terms of the overall computational cost. Ourown experience is that (GO–MARK4) is a competitive strategy in every example that hasbeen tested. Consequently, we have made it the default option within the code. In order to demonstrate the effectiveness of the goal-oriented adaptive strategy describedin the previous subsection, let us consider the following test example.
Example 3.
Let us consider the model problem given by (2.4) with A = Id on a slitdomain D δ = ( − , \ T δ , where T δ = conv { (0 , , ( − , δ ) , ( − , − δ ) } with δ = 0 . u to the (primal) problem in this example exhibits a singularityinduced by the slit in the domain (see Figure 7(b)). Our aim, however, is to demonstratethe capability of the software to approximate the value of u at some fixed point x ∈ D away from the slit (in the experiments below, we set x = (0 . , − . − (a) − − . (b)Figure 7. Example 3: (a) initial coarse mesh in the GOAFEM algorithm; (b) the primal Galerkin solution. the corresponding bounded goal functional G , it is common to fix a sufficiently small r > mollifier g as follows (cf. [48]): g ( x ) = g ( x ; x , r ) := (cid:40) C exp (cid:16) − r r −| x − x | (cid:17) if | x − x | < r, C is chosen such that (cid:82) D g ( x ) d x = 1 (for sufficiently small r suchthat supp( g ( x ; x , r )) ⊂ D , one has C ≈ . r − ; see, e.g., [48]). Then, the functional G in (3.2) reads as G ( v ) = (cid:90) D g ( x ) v ( x ) d x for all v ∈ H ( D ) . Note that if u ( x ) is continuous in a neighborhood of x , then G ( u ) converges to the pointvalue u ( x ) as r tends to zero.We started with the coarse triangulation depicted in Figure 7(a) in all the experi-ments. In the first experiment, we fixed the stopping tolerance to be tol = 3e-4 and ran r = 0 . , T = 5653 r = 0 . , T = 3834 r = 0 . , T = 2815 − − . . . . − − . . . . − − . . . . Figure 8. Example 3: adaptively refined triangulations (top row) and the dual Galerkin solutions (bottomrow) computed using the mollifier g in (3.5) with r = 0 . , . , . − − − − − O ( N − . ) O ( N − )degrees of freedom e rr o r e s t i m a t e r = 0 . | G ( u ref ) − G ( u h ) | (a) . . . e ff ec t i v i t y i nd e x r = 0 . (b)Figure 9. Example 3: (a) error estimates at each iteration of the GOAFEM algorithm employing themarking strategy (GO–MARK4); (b) the associated effectivity indices. the GOAFEM algorithm to compute dual Galerkin solutions for different values of theradius r in (3.5). For the SOLVE step, we used P approximations for both primal anddual solutions. Within the ESTIMATE module, the energy errors in both solutions wereestimated using the two-level error estimation strategy (EES3) described earlier. Giventhe error indicators for primal and dual solutions, the algorithm employed the edge-basedD¨orfler marking (2.3) with θ = 0 . r = 0 . , . , .
5. We note that the triangulations generatedby the algorithm adapt to the features of both primal and dual solutions: the triangula-tions are refined in the vicinity of each corner (with particularly strong refinement nearthe tip of the slit) and in a neighborhood of x (with stronger refinement for smallervalues of r ).Focusing now on the case r = 0 .
2, we set tol = 8e-5 in the second experiment andran the GOAFEM algorithm without changing the settings. The results we obtained areshown in Figure 9. In Figure 9(a), we plot the energy error estimates for primal anddual Galerkin approximations, the estimates of the error in the goal functional, as well asthe reference errors in the goal functional (i.e., | G ( u ref ) − G ( u h ) | ) at each iteration of theGOAFEM algorithm. Here, the reference Galerkin solution u ref is computed using thetriangulation obtained by two uniform refinements of the final triangulation generated bythe GOAFEM algorithm. We observe that all error estimates as well as the reference errorin the goal functional converge with optimal rates. The effectivity indices for the goal-oriented error estimation at each iteration of the algorithm are plotted in Figure 9(b). Thisplot shows that the product of energy error estimates for the primal and dual Galerkinsolutions provides a reasonably accurate estimate for the error in approximating the goalfunctional G ( u ). In this section, we address the design of components of the adaptive algorithm whenworking in a parametric PDE setting. Before discussing the details of our implementation,13e formulate the model problem with parametric input data and briefly describe theidea of stochastic Galerkin FEM (SGFEM). Readers interested in theoretical aspects ofSGFEM are referred to [31], [19] and [2].
Let D ⊂ R be a Lipschitz domain (called the physical domain) with polygonal boundary ∂D and let Γ := (cid:81) ∞ m =1 [ − ,
1] denote the infinitely-dimensional hypercube (called theparameter domain). We consider the elliptic boundary value problem −∇ · ( a ∇ u ) = f in D × Γ ,u = 0 on ∂D × Γ , (4.1)where the scalar coefficient a (and, hence, the solution u ) depends on a countably infinitenumber of scalar parameters, i.e., a = a ( x , y ) and u = u ( x , y ) with x ∈ D , y ∈ Γ, and thedifferentiation in ∇ is with respect to x = ( x , x ). We assume that f = f ( x ) ∈ H − ( D )and that the parametric coefficient a has affine dependence on the parameters, i.e., a ( x , y ) = a ( x ) + ∞ (cid:88) m =1 y m a m ( x ) for x ∈ D and y = ( y m ) m ∈ N ∈ Γ , (4.2)where the scalar functions a m ∈ W , ∞ ( D ) ( m ∈ N ) satisfy the following inequalities:0 < a min0 ≤ a ( x ) ≤ a max0 < ∞ for almost all x ∈ D (4.3)and τ := 1 a min0 ∞ (cid:88) m =1 (cid:107) a m (cid:107) L ∞ ( D ) < . (4.4)The weak formulation of (4.1) is posed in the framework of the Bochner space V := L π (Γ; H ( D )). Here, π = π ( y ) is a probability measure on (Γ , B (Γ)) with B (Γ) beingthe Borel σ -algebra on Γ, and we assume that π ( y ) is the product of symmetric Borelprobability measures π m on [ − , π ( y ) = (cid:81) ∞ m =1 π m ( y m ). For a given f ∈ H − ( D ),the weak solution u ∈ V satisfies B ( u, v ) = F ( v ) := (cid:90) Γ (cid:90) D f ( x ) v ( x , y ) d x d π ( y ) for all v ∈ V. (4.5)Here, B ( u, v ) := B ( u, v ) + ∞ (cid:88) m =1 (cid:90) Γ (cid:90) D y m a m ( x ) ∇ u ( x , y ) · ∇ v ( x , y ) d x d π ( y ) ,B ( u, v ) := (cid:90) Γ (cid:90) D a ( x ) ∇ u ( x , y ) · ∇ v ( x , y ) d x d π ( y ) . (4.6)The following three coefficient expansions (CEs) of the type (4.2) are implemented inStochastic T-IFISS.(CE1) Let D = ( − , and suppose that the coefficient a = a ( x , y ) in (4.1) is aparametric representation of a second-order random field with prescribed mean E [ a ] and14ovariance function Cov[ a ]. Assume that Cov[ a ] is the separable exponential covariancefunction given by Cov[ a ]( x , x (cid:48) ) = σ exp (cid:18) − | x − x (cid:48) | l − | x − x (cid:48) | l (cid:19) , where x , x (cid:48) ∈ D , σ denotes the standard deviation, and l , l are correlation lengths.In this case, a ( x , y ) can be written in the form (4.2) using the Karhunen–Lo`eve-typeexpansion [40] such that a ( x ) := E [ a ]( x ) , a m ( x ) := c (cid:112) λ m ϕ m ( x ) , m ∈ N , where { ( λ m , ϕ m ) } ∞ m =1 are the eigenpairs of the operator (cid:82) D Cov[ a ]( x , x (cid:48) ) ϕ ( x (cid:48) )d x (cid:48) , and theconstant c > c y m ) = 1 for all m ∈ N . Note that analyticalexpressions for λ m and ϕ m for the square domain D follow from the corresponding formulasderived in [31, pp. 28–29] for the one-dimensional case.(CE2) Following [21, Section 11.1], we set a ( x ) := 1, x ∈ D , and choose the coeffi-cients a m ( x ) in (4.2) to represent planar Fourier modes of increasing total order: a m ( x ) := α m cos(2 πβ ( m ) x ) cos(2 πβ ( m ) x ) for all m ∈ N , (4.7)where α m := Am − ˜ σ are the amplitudes of the coefficients, ˜ σ >
1, the constant A is chosensuch that τ = Aζ (˜ σ ) = 0 . ζ denotes the Riemann zeta function), and β , β aredefined as β ( m ) := m − k ( m )( k ( m ) + 1) / β ( m ) := k ( m ) − β ( m ) , with k ( m ) := (cid:98)− / (cid:112) / m (cid:99) for all m ∈ N . The following two values of the decayparameter ˜ σ are used in test problems: ˜ σ = 2 (slow decay) and ˜ σ = 4 (fast decay).(CE3) Finally, we consider the following parametric coefficient (cf. [40, Example 9.37]): a ( x , ˜ y ) = 1 + c ∞ (cid:88) i =0 ∞ (cid:88) j =0 (cid:112) λ ij ϕ ij ( x ) y ij , ˜ y = ( y ij ) i,j ∈ N , (4.8)where λ ij = ¯ λ i ¯ λ j , ϕ ij ( x ) = ¯ ϕ i ( x ) ¯ ϕ j ( x ), y ij ∈ [ − ,
1] ( i, j ∈ N ) with ¯ λ := 1 / ϕ ( t ) := 1, ¯ λ k := exp( − πk (cid:96) ), ¯ ϕ k ( t ) := √ πkt ) ( k ∈ N ), and the constant c > c y i,j ) = 1 for all i, j ∈ N .As shown in [40, Example 9.37], the parametric representation (4.8) stems from theKarhunen-Lo`eve expansion of a random field with the mean E [ a ] = 1 and covariancefunction close to the isotropic covariance c ( x ) = (4 (cid:96) ) − exp( − π ( x + x ) / (4 (cid:96) )), where (cid:96) is the correlation length. For our implementation, we set (cid:96) = 1 and reorder the terms inthe double sum in (4.8) to write the coefficient a ( x , ˜ y ) in the form (4.2) with a ( x ) := 1 , a m ( x ) := c (cid:112) λ m ϕ m ( x ) , y m ∈ [ − , , m ∈ N , where λ ≥ λ ≥ . . . .In each coefficient expansion (CE1)–(CE3), parameters y m ( m ∈ N ) are the imagesof independent mean-zero random variables on Γ m = [ − , π m = d y m / L π m (Γ m ) is comprised of scaled Legendre polynomials.15RV2) Truncated Gaussian random variables. In this case, d π m = p ( y m )d y m with p ( y m ) = exp( − y m / (2 σ )) σ √ π (2Φ(1 /σ ) − , m ∈ N , (4.9)where Φ( · ) is the Gaussian cumulative distribution function and σ is a parameter measur-ing the standard deviation. The corresponding orthonormal polynomial basis in L π m (Γ m )is formed by the so-called Rys polynomials; see, e.g., [30, Example 1.11].A key observation that motivates the stochastic Galerkin FEM is that the Bochnerspace V = L π (Γ; H ( D )) is isometrically isomorphic to H ( D ) ⊗ L π (Γ). Mimicking thistensor-product construction, the finite-dimensional subspace V XP ⊂ V is defined as V XP := X ⊗ P ; here, X = X h is a finite element space associated with a conforming triangulation T h of D and P = P P is a polynomial space on Γ associated with a finite index set P .Specifically, X = X h := span { φ i ; i = 1 , . . . N X } ⊂ H ( D ) , N X := dim( X ) (4.10)and P = P P := span (cid:110) P ν ( y ) = (cid:89) m ∈ N P mν m ( y m ); ν ∈ P (cid:111) ⊂ L π (Γ) with P ⊂ I , where { P mn : n ∈ N } is an orthonormal basis of L π m ( − ,
1) and I denotes the countableset of finitely supported multi-indices, i.e., I := (cid:8) ν = ( ν m ) m ∈ N ; ν m ∈ N for all m ∈ N , ν ) < ∞ (cid:9) with supp( ν ) := { m ∈ N ; ν m (cid:54) = 0 } .The Galerkin discretization of (4.5) reads as follows: find u XP ∈ V XP such that B ( u XP , v ) = F ( v ) for all v ∈ V XP . (4.11)Note that dim( V XP ) = dim( X h ) × dim( P P ). Therefore, if a large number of random vari-ables is used to represent the input data, then computing high-fidelity stochastic Galerkinapproximations with standard polynomial subspaces on Γ (e.g., the spaces of tensor-product or complete polynomials) becomes prohibitively expensive. This motivates thedevelopment of adaptive SGFEM algorithms that incrementally refine spatial ( X -) andparametric ( P -) components of Galerkin approximations by iterating the standard adap-tive loop (2.1). The implementation details are discussed in the following subsections. SOLVE . Linear algebra aspects of SGFEM
Recalling the definitions of X h and P P , the Galerkin solution u XP is sought in the form u XP = N X (cid:88) i =1 N P (cid:88) j =1 u ij φ i ( x ) P κ ( j ) ( y ) , (4.12)where N P := dim( P ) = P , κ is a bijection { , , . . . , N P } → P , and the coefficients u ij are computed by solving the linear system A u = b with block structure. Specifically, thesolution vector u and the right-hand side vector b are given by u = [ u u . . . u N P ] T and b = [ b b . . . b N P ] T , The module
SOLVE is implemented for both P and P finite element approximations in the currentversion of the software, whereas the error estimation module (and hence, the adaptive algorithm) is onlyimplemented for P approximation. u j := [ u j u j . . . u N X j ] T , j = 1 , . . . , N P , [ b t ] s := (cid:104) , P κ ( t ) (cid:105) π (cid:90) D f ( x ) φ s ( x ) d x , s = 1 , . . . , N X , t = 1 , . . . , N P ;the coefficient matrix A is given by (see, e.g., [40, Section 9.5]) A = G ⊗ K + M P (cid:88) m =1 G m ⊗ K m , (4.13)where M P is the number of active parameters in the index set P ,[ G ] tj := (cid:104) P κ ( j ) , P κ ( t ) (cid:105) π = δ tj , [ G m ] tj := (cid:104) y m P κ ( j ) , P κ ( t ) (cid:105) π ( m = 1 , . . . , M P )with t, j = 1 , . . . , N P , and K m are the finite element (stiffness) matrices defined by[ K m ] si := (cid:90) D a m ∇ φ i · ∇ φ s d x , m = 0 , , . . . , M P , s, i = 1 , . . . , N X . The design of an efficient linear solver is a crucial ingredient of the stochastic Galerkinapproximation process. Rather than computing a (memory intensive) sparse factorizationof the coefficient matrix, a matrix-free iterative solver is needed. The key idea is that thematrix-vector products with A can be computed from its sparse matrix components byexploiting the Kronecker product structure, without assembling A itself. The iterativesolver that enables this process within T-IFISS is a bespoke implementation of the Mini-mum Residual algorithm, called EST MINRES [53]. The MINRES algorithm is designedto solve symmetric (possibly indefinite) linear equation systems and requires the actionof A on a given vector at every iteration, see [25, Section 2.4]. Using this strategy thestorage overhead (in addition to the component matrices K , . . . , K M P , G , . . . , G M P ) isfor five vectors of length N P · N X .A crucial ingredient in the design of a fast iterative solver is preconditioning . Thestandard choice of preconditioning operator in this context is the parameter-free matrixoperator P = G ⊗ K = I ⊗ K . The action of P − r , where r is the current residual vector, is needed at every iteration—this can be done efficiently by computing a single sparse triangular factorization of thematrix K and then performing N P forward and backward substitutions on the compo-nents of the residual vector. Theoretical analysis of the preconditioned operator givenin [47] shows that the eigenvalues of the preconditioned operator are bounded away fromzero and bounded away from infinity independently of the discretization parameters N X and N P . This means that the number of preconditioned EST MINRES iterations neededto satisfy a fixed residual reduction tolerance will not grow unboundedly when the dis-cretization parameters are changed. In practice, the number of iterations needed to satisfythe default tolerance of is less than 20, independent of the finite element mesh res-olution and the number of active indices. ESTIMATE . Error estimation in SGFEM
Stochastic Galerkin approximations are built from two distinct discretizations: the spatial(finite element) discretization over the physical domain D and the parametric (polyno-mial) approximation on the parameter domain Γ. Therefore, there are two distinct sources17f discretization error arising from the choice of the finite element space X and the poly-nomial space P . This fact determines the structure of a posteriori estimates for the energyerrors in SGFEM approximations as combinations of spatial and parametric contributions(cf. [21, 22, 10, 15]).The spatial errors in SGFEM approximations are estimated by extending the strategies(EES1)–(EES3) described in § e X | K ( K ∈ T h ), now lives in thetensor-product space Y | K ⊗ P P , where Y | K is the local space of piecewise linear or piecewisequadratic bubble functions. These error estimators are computed by solving local residualproblems of the following type (see [10, Section 6.2] for details): B ,K ( e X | K , v ) = Res K ( a, f, u XP ; v ) ∀ v ∈ Y | K ⊗ P P , (4.14)where B ,K is the elementwise bilinear form associated with the parameter-free term a in the coefficient expansion (cf. (4.6)). This construction of the error estimator enablesfast linear algebra for solving (4.14). Indeed, the coefficient matrix in the linear systemassociated with (4.14) has a very simple structure: it is the Kronecker product of a3 × N P = dim( P ). As aresult, the action of the inverse of this coefficient matrix can be effected by a block LDL T factorization of the element stiffness matrices followed by a sequence of N P backwardand forward substitutions. Furthermore, since the factorizations and triangular solvesare logically independent, the entire computation is vectorized over the finite elementsthat define the spatial subdivision. We refer to [12, Section 3] for details of the globalhierarchical (EES2) and the two-level (EES3) error estimation strategies in the contextof the SGFEM.The parametric errors in SGFEM approximations are estimated using the hierarchicalapproach in the spirit of [6]. To that end, we first introduce the finite index set Q P as a“neighborhood” of the index set P . More precisely, for a fixed M ∈ N , we define Q P := (cid:8) ν ∈ I \ P ; ν = µ ± ε ( m ) for some µ ∈ P and some m = 1 , . . . , M P + M (cid:9) , (4.15)where ε ( m ) := ( ε ( m )1 , ε ( m )2 , . . . ) ( m ∈ N ) denotes the Kronecker delta index such that ε ( m ) k = δ mk for all k ∈ N , and M P ∈ N is the number of active parameters in P .For a given P ⊂ I , the index set Q P contains only those “neighbors” of all indices in P that have up to M P + M active parameters, that is M parameters more than currentlyactivated in the index set P (we refer to [15, Section 4.2] for theoretical underpinnings ofthis construction). Then, the parametric error estimator e P is computed as a combinationof the contributing estimators e ( ν ) P associated with individual indices ν ∈ Q P , i.e., e P = (cid:80) ν ∈ Q P e ( ν ) P , where each contributing estimator e ( ν ) P ∈ X ⊗ span( P ν ), ν ∈ Q P , is computedby solving the linear system associated with the following discrete formulation: B ( e ( ν ) P , v P ν ) = F ( v P ν ) − B ( u XP , v P ν ) for all v ∈ X . (4.16)Note that the coefficient matrix of this linear system represents the assembled stiffnessmatrix corresponding to the parameter-free term a in (4.2), and is therefore the samefor all ν ∈ Q P . Once the stiffness matrix has been factorized, the estimators e ( ν ) P arecomputed independently by using forward and backward substitutions.Once the spatial and parametric error estimators have been computed, the total errorestimate η is calculated via η := (cid:16) (cid:107) e X (cid:107) + (cid:107) e P (cid:107) (cid:17) / = (cid:18) (cid:88) K ∈T h (cid:13)(cid:13) e X | K (cid:107) ,K + (cid:88) ν ∈ Q P (cid:13)(cid:13) e ( ν ) P (cid:13)(cid:13) (cid:19) / , (4.17)18here (cid:107) · (cid:107) (resp., (cid:107) · (cid:107) ,K ) denotes the norm induced by the bilinear form B (resp., B ,K ). The module
ESTIMATE supplies local spatial error indicators associated with elements oredges of triangulation (e.g., e X | K for the error estimation strategy (EES1)) as well as thecontributing parametric error indicators e ( ν ) P associated with individual indices ν ∈ Q P .In the module MARK , the largest error indicators are selected independently for spatialand for parametric components of Galerkin approximations. To that end, one of themarking strategies described in § θ X and θ P , respectively. However,a simple modification of the code will allow one to use different marking strategies fordifferent components of Galerkin approximations.Thus, at each iteration of the adaptive SGFEM algorithm, the output of the module MARK contains two sets: the set of marked elements in the current mesh T h to be refined(or, the set of edges to be bisected) and the set M ⊆ Q P of marked indices to be addedto the current index set P (note that choosing M > V XP is then enhanced within the module REFINE by performing eitherspatial refinement (as described in §
2) or parametric refinement (simply by adding M to P ). The question then arises which type of refinement (spatial or parametric) should beperformed at a given iteration.A traditional strategy for choosing between the two refinements is based on thedominant error estimator contributing to the total error estimate η defined in (4.17);cf. [21, 22, 15]. This strategy is referred to as version 1 of the adaptive algorithm imple-mented in Stochastic T-IFISS. An alternative strategy is referred to as version 2 of theimplemented algorithm: here, the refinement type that leads to a larger estimated errorreduction is chosen at each iteration; see [13, 11]. This strategy exploits the fact that local spatial error indicators (e.g., (cid:107) e X | K (cid:107) ,K ( K ∈ T h ) in the error estimation strategy (EES1))and individual parametric error indicators (cid:107) e ( ν ) P (cid:107) ( ν ∈ Q P ) provide effective estimates ofthe error reduction that would be achieved by performing, respectively, a local refine-ment of the current mesh (e.g., by refining the element K ) and a selective enrichment ofthe parametric component of the current Galerkin approximation (by adding the index ν ∈ Q P to the current index set P ). We refer to [10, Theorem 5.1] and [11, Corollary 3]for the underpinning theoretical results and to [13] and [11, Section 5] for comprehensivenumerical studies of the two versions of the adaptive algorithm and different markingstrategies. We conclude this section with a representative case study that demonstrates the efficiencyof our adaptive SGFEM algorithm.
Example 4.
We consider the parametric model problem (4.1) on the L-shaped domain D = ( − , \ ( − , . We set f ( x ) = (1 − x ) − . and choose the parametric coefficient a ( x , y ) in the form (4.2), where a and a m ( m ∈ N ) are as specified by the coefficientexpansion (CE2) with ˜ σ = 2 and y m ( m ∈ N ) are the images of independent truncatedGaussian random variables with zero mean on Γ m = [ − ,
1] (see (RV2), where we set19 − . . (a) − − · − (b) − − (c)Figure 10. Example 4: (a)–(b) mean and variance of the SGFEM solution; (c) a typical locally refinedmesh generated by the adaptive SGFEM algorithm. σ = 1). The mean and the variance of an SGFEM solution to this problem are shown inFigure 10(a)–(b), whereas Figure 10(c) depicts a typical locally refined mesh generatedby the adaptive SGFEM algorithm. Note how the adaptively refined mesh effectivelyidentifies the areas of singular behavior of the mean field—in the vicinity of the reentrantcorner (due to a geometric singularity) and in the vicinity of the edge x = 1 (due to asingular right-hand side function).When running the adaptive SGFEM algorithm, we start with a uniform mesh consist-ing of 96 right-angled triangles and an initial index set P := { (0 , , , . . . , ) , (1 , , , . . . , ) } .For the error estimation module, we employ the two-level spatial error estimator (EES3)combined with the hierarchical parametric error estimators associated with individualindices ν ∈ Q P (see (4.16) and (4.15), where we set M := 1). We use D¨orfler mark-ing for spatial and parametric components of Galerkin approximations (for the spatialcomponent, we mark the edges of the mesh).In our first experiment we ran the adaptive SGFEM algorithm (in version 2) with thefollowing four sets of D¨orfler marking parameters (with different stopping tolerances):(i) θ X = 1, θ P = 1 (no adaptivity in either of the components), tol=5.1e-3 ;(ii) θ X = 0 . θ P = 1 (adaptive refinement only for spatial component), tol=4e-3 ; − − − degrees of freedom e rr o r e s t i m a t e θ X = 1, θ P = 1 θ X = 0 . θ P = 1 θ X = 1, θ P = 0 . θ X = 0 . θ P = 0 . O ( N − / ) Figure 11. Example 4: error estimates at each iteration of the adaptive SGFEM algorithm for differentsets of marking parameters.
Table 3. Example 4: evolution of the index set when running adaptive SGFEM algorithm with θ X = 0 . θ P = 0 . (iii) θ X = 1, θ P = 0 . tol=4e-3 ;(iv) θ X = 0 . θ P = 0 . tol=3e-3 .For each run of the algorithm, the error estimates computed at each iteration are plottedin Figure 11. The error estimates can be seen to decrease at every iteration. However,in cases (i)–(iii), the decay rate either eventually deteriorates (cases (i) and (ii)), due tothe number of degrees of freedom growing very fast, or it is significantly slower (case (iii))than in the case of adaptivity being used for both components of SGFEM approximations(case (iv)). This shows that, for the same level of accuracy, adaptive refinement of both components results in more balanced approximations with fewer degrees of freedom andleads to the fastest convergence rate for the parameter choices in this experiment.To give an indication of the algorithmic efficiency, we provide further details of the runin case (iv). In this computation, which took 927 seconds, the stopping tolerance was met after 19 iterations. The final triangulation generated by the algorithm comprised545,636 finite elements with 271,599 interior vertices (the latter number defines the dimen-sion of the corresponding finite element space). For the final polynomial approximationon Γ, the algorithm produced an index set P of cardinality 23 with 6 active parameters;the evolution of the index set throughout the computation is shown in Table 3. The totalnumber of d.o.f. in the SGFEM solution at the final iteration was equal to 6,246,777. In21 − − − degrees of freedom e rr o r e s t i m a t e Adaptive SGFEM ( θ X = 0 . θ P = 0 . η (total) (cid:107) e X (cid:107) (spatial) (cid:107) e P (cid:107) (parametric) (a) . . . . degrees of freedom e ff ec t i v i t y i nd e x Adaptive SGFEM ( θ X = 0 . θ P = 0 . Figure 12(a), we show the interplay between the spatial and parametric contributions tothe total error estimates η (see (4.17)) at each iteration of the adaptive algorithm. Theeffectivity indices for the total error estimates are plotted in Figure 12(b). These werecalculated using the reference Galerkin solution computed with P approximations onthe final triangulation generated by the algorithm and with the polynomial space P P ∪ Q P ,where P is the final index set produced by the algorithm and Q P is the “neighborhood”of P as defined in (4.15) with M = 1 (the total number of d.o.f. in this reference solutionwas 37,020,322). The T-IFISS software framework provides many opportunities for experimentation andexploration. It is also an invaluable teaching tool for numerical analysis and computationalengineering courses with an emphasis on contemporary finite element analysis. StochasticT-IFISS has been recently extended to incorporate the goal-oriented error estimation andadaptivity in the context of stochastic Galerkin approximations for parametric ellipticPDEs (see [12] for details of the algorithm and numerical results). The toolbox hasalso been used for numerical testing of the adaptive algorithm proposed for parameter-dependent linear elasticity problems; see [37, 36]. Future developments would include theextension to problems with non-affine parametric representations of inputs (see [16]) andthe implementation of the multilevel adaptive SGFEM algorithm in the spirit of [21, 18].
Acknowledgement.
The authors are grateful to Qifeng Liao (ShanghaiTech University)for his inputs to T-IFISS project at its early stages. The authors would also like tothank Dirk Praetorius and Michele Ruggeri (both at Technical University of Vienna) fortheir contributions to the development of the toolbox components for goal-oriented errorestimation and adaptivity. 22 eferences [1]
M. Ainsworth and J. T. Oden , A posteriori error estimation in finite elementanalysis , Pure and Applied Mathematics (New York), Wiley, 2000.[2]
I. Babuˇska, R. Tempone, and E. Zouraris , Galerkin finite element approxima-tions of stochastic elliptic partial differential equations , SIAM J. Numer. Anal., 42(2004), pp. 800–825.[3]
I. Babuˇska and M. Vogelius , Feedback and adaptive finite element solution ofone-dimensional boundary value problems. , Numer. Math., 44 (1984), pp. 75–102.[4]
W. Bangerth, R. Hartmann, and G. Kanschat , deal.II – a general purpose ob-ject oriented finite element library , ACM Trans. Math. Software, 33 (2007), pp. 24/1–24/27. .[5] R. E. Bank , PLTMG: a software package for solving elliptic partial differentialequations. Users’ guide 8.0 , vol. 5 of Software, Environments, and Tools, Societyfor Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1998. .[6]
R. E. Bank and A. Weiser , Some a posteriori error estimators for elliptic partialdifferential equations , Math. Comp., 44 (1985).[7]
E. B¨ansch , Local mesh refinement in 2 and 3 dimensions , IMPACT Comput. Sci.Engrg., 3 (1991), pp. 181–191.[8]
R. Becker, E. Estecahandy, and D. Trujillo , Weighted marking for goal-oriented adaptive finite element methods , SIAM J. Numer. Anal., 49 (2011), pp. 2451–2469.[9]
R. Becker and R. Rannacher , An optimal control approach to a posteriori errorestimation in finite element methods , Acta Numer., 10 (2001), pp. 1–102.[10]
A. Bespalov, C. E. Powell, and D. Silvester , Energy norm a posteriorierror estimation for parametric operator equations , SIAM J. Sci. Comput., 36 (2014),pp. A339–A363.[11]
A. Bespalov, D. Praetorius, L. Rocchi, and M. Ruggeri , Convergence ofadaptive stochastic Galerkin FEM , SIAM J. Numer. Anal., 57 (2019), pp. 2359–2382.[12]
A. Bespalov, D. Praetorius, L. Rocchi, and M. Ruggeri , Goal-orientederror estimation and adaptivity for elliptic PDEs with parametric or uncertain inputs ,Comput. Methods Appl. Mech. Engrg., 345 (2019), pp. 951–982.[13]
A. Bespalov and L. Rocchi , Efficient adaptive algorithms for elliptic PDEs withrandom data , SIAM/ASA J. Uncertain. Quantif., 6 (2018), pp. 243–272.[14] ,
Stochastic T-IFISS , February 2019. Available online at http://web.mat.bham.ac.uk/A.Bespalov/software/index.html .[15]
A. Bespalov and D. Silvester , Efficient adaptive stochastic Galerkin methods forparametric operator equations , SIAM J. Sci. Comput., 38 (2016), pp. A2118–A2140.2316]
A. Bespalov and F. Xu , A posteriori error estimation and adaptivity in stochas-tic Galerkin FEM for parametric elliptic PDEs: beyond the affine case . Preprint,arXiv:1903.06520 [math.NA], 2019.[17]
M. Blatt, A. Burchardt, A. Dedner, C. Engwer, J. Fahlke, B. Fle-misch, C. Gersbacher, C. Gr¨aser, F. Gruber, C. Gr¨uninger, D. Kempf,R. Kl¨ofkorn, T. Malkmus, S. M¨uthing, M. Nolte, M. Piatkowski, andO. Sander , The distributed and unified numerics environment, version 2.4 , Arch.Num. Soft., 4 (2016), pp. 13–29. .[18]
A. J. Crowder, C. E. Powell, and A. Bespalov , Efficient adaptive multilevelstochastic Galerkin approximation using implicit a posteriori error estimation , SIAMJ. Sci. Comput., 41 (2019), pp. A1681–A1705.[19]
M. K. Deb, I. Babuˇska, and J. T. Oden , Solution of stochastic partial differ-ential equations using Galerkin finite element techniques , Comput. Methods Appl.Mech. Engrg., 190 (2001), pp. 6359–6372.[20]
W. D¨orfler , A convergent adaptive algorithm for Poisson’s equation , SIAM J.Numer. Anal., 33 (1996), pp. 1106–1124.[21]
M. Eigel, C. J. Gittelson, C. Schwab, and E. Zander , Adaptive stochasticGalerkin FEM , Comput. Methods Appl. Mech. Engrg., 270 (2014), pp. 247–269.[22] ,
A convergent adaptive stochastic Galerkin finite element method with quasi-optimal spatial meshes , ESAIM Math. Model. Numer. Anal., 49 (2015), pp. 1367–1398.[23]
M. Eigel and E. Zander , ALEA – A python framework for spectral methods andlow-rank approximations in uncertainty quantification . https://bitbucket.org/aleadev/alea/src .[24] H. Elman, A. Ramage, and D. Silvester , IFISS: a computational laboratory forinvestigating incompressible flow problems , SIAM Review, 56 (2014), pp. 261–273.[25]
H. Elman, D. Silvester, and A. Wathen , Finite Elements and Fast Itera-tive Solvers: with Applications in Incompressible Fluid Dynamics , Oxford UniversityPress, Oxford, UK, 2014. Second Edition, xiv+400 pp. ISBN: 978-0-19-967880-8.[26]
H. C. Elman, A. Ramage, and D. J. Silvester , Algorithm 886: IFISS, aMatlab toolbox for modelling incompressible flow , ACM Trans. Math. Software, 33(2007), article 14.[27]
K. Eriksson, D. Estep, P. Hansbo, and C. Johnson , Introduction to adaptivemethods for differential equations , Acta Numer., 4 (1995), pp. 105–158.[28]
M. Feischl, D. Praetorius, and K. G. van der Zee , An abstract analysis ofoptimal goal-oriented adaptivity , SIAM J. Numer. Anal., 54 (2016), pp. 1423–1448.[29]
S. Funken, D. Praetorius, and P. Wissgott , Efficient implementation ofadaptive P1-FEM in Matlab , Comput. Methods Appl. Math., 11 (2011), pp. 460–490. .2430]
W. Gautschi , Orthogonal polynomials: computation and approximation , NumericalMathematics and Scientific Computation, Oxford University Press, New York, 2004.[31]
R. G. Ghanem and P. D. Spanos , Stochastic finite elements: a spectral approach ,Springer-Verlag, New York, 1991.[32]
M. B. Giles and E. S¨uli , Adjoint methods for PDEs: a posteriori error analysisand postprocessing by duality , Acta Numer., 11 (2002), pp. 145–236.[33]
A. Gopal and L. N. Trefethen , New Laplace and Helmholtz solvers , Proc. Nat.Acad. Sci., 116 (2019), pp. 10223–10225.[34]
F. Hecht , New development in FreeFem++ , J. Numer. Math., 20 (2012), pp. 251–265. https://freefem.org/ .[35]
M. Holst and S. Pollock , Convergence of goal-oriented adaptive finite elementmethods for nonsymmetric problems , Numer. Methods Partial Differential Equations,32 (2016), pp. 479–509.[36]
A. Khan, A. Bespalov, C. E. Powell, and D. J. Silvester , Robust a pos-teriori error estimation for stochastic Galerkin formulations of parameter-dependentlinear elasticity equations . Preprint, arXiv:1810.07440 [math.NA], 2018.[37]
A. Khan, C. E. Powell, and D. J. Silvester , Robust preconditioning forstochastic Galerkin formulations of parameter-dependent nearly incompressible elas-ticity equations , SIAM J. Sci. Comput., 41 (2019), pp. A402–A421.[38]
I. Kossaczky , A recursive approach to local mesh refinement in two and three di-mensions , J. Comput. Appl. Math., 55 (1995), pp. 275–288.[39]
A. Logg, K.-A. Mardal, G. N. Wells, et al. , Automated Solution ofDifferential Equations by the Finite Element Method , Springer, 2012. https://fenicsproject.org/ .[40]
G. J. Lord, C. E. Powell, and T. Shardlow , An introduction to computationalstochastic PDEs , Cambridge Texts in Applied Mathematics, Cambridge UniversityPress, New York, 2014.[41]
W. F. Mitchell , A comparison of adaptive refinement techniques for elliptic prob-lems , ACM Trans. Math. Software, 15 (1989), pp. 326–347.[42]
M. S. Mommer and R. Stevenson , A goal-oriented adaptive finite elementmethod with convergence rates , SIAM J. Numer. Anal., 47 (2009), pp. 861–886.[43]
P. Mund and E. P. Stephan , An adaptive two-level method for the coupling ofnonlinear FEM-BEM equations , SIAM J. Numer. Anal., 36 (1999), pp. 1001–1021.[44]
P. Mund, E. P. Stephan, and J. Weiße , Two-level methods for the single layerpotential in R , Computing, 60 (1998), pp. 243–266.[45] R. H. Nochetto and A. Veeser , Primer of adaptive finite element methods ,in Multiscale and Adaptivity: Modeling, Numerics and Applications, vol. 2040,Springer-Verlag Berlin Heidelberg, 2012, pp. 125–225.2546]
P.-O. Persson and G. Strang , A simple mesh generator in Matlab , SIAM Rev.,46 (2004), pp. 329–345. http://persson.berkeley.edu/distmesh/ .[47]
C. E. Powell and H. C. Elman , Block-diagonal preconditioning for spectralstochastic finite-element systems , IMA Journal of Numerical Analysis, 29 (2009),pp. 350–375.[48]
S. Prudhomme and J. T. Oden , On goal-oriented error estimation for ellipticproblems: application to the control of pointwise errors , Comput. Methods Appl.Mech. Engrg., 176 (1999), pp. 313–331. New advances in computational methods(Cachan, 1997).[49]
M. C. Rivara , Mesh refinement processes based on the generalized bisection of sim-plices , SIAM J. Numer. Anal., 21 (1984), pp. 604–613.[50]
L. Rocchi , Adaptive algorithms for partial differential equations with parametricuncertainty , PhD thesis, University of Birmingham, 2019. Electronically publishedat https://etheses.bham.ac.uk/id/eprint/9157/.[51]
A. Schmidt and K. G. Siebert , Design of adaptive finite element software. Thefinite element toolbox ALBERTA , vol. 42 of Lecture Notes in Computational Scienceand Engineering, Springer-Verlag, Berlin, 2005. .[52]
E. G. Sewell , Automatic generation of triangulations for piecewise polynomial ap-proximation , PhD thesis, Purdue University, 1972.[53]
D. J. Silvester and V. Simoncini , An optimal iterative solver for symmetricindefinite systems stemming from mixed approximation , ACM Trans. Math. Software,37 (2011), pp. 42/1–42/22.[54]
R. Stevenson , The completion of locally refined simplicial partitions created bybisection , Math. Comp., 77 (2008), pp. 227–241.[55]
L. N. Trefethen , Posting on NA Digest at (29 November 2018).[56]
E. Zander , SGLib v0.9 . https://github.com/ezander/sglibhttps://github.com/ezander/sglib