Automatically Stable Discontinuous Petrov-Galerkin Methods for Stationary Transport Problems: Quasi-Optimal Test Space Norm Antti H. Niemi, Nathaniel O. Collier, Victor M. Calo ∗ † ‡ King Abdullah University of Science and Technology (KAUST) 2 Thuwal 23955-6900, Kingdom of Saudi Arabia 1 0 2 n a Abstract J 9 We investigate the application of the discontinuous Petrov-Galerkin (DPG) finite el- ement framework to stationary convection-diffusion problems. In particular, we demon- ] A strate how the quasi-optimal test space norm can be utilized to improve the robustness N of the DPG method with respect to vanishing diffusion. We numerically compare coarse- . mesh accuracy of the approximation when using the quasi-optimal norm, the standard h norm, and the weighted norm. Our results show that the quasi-optimal norm leads to t a more accurate results on three benchmark problems in two spatial dimensions. We ad- m dress the problems associated to the resolution of the optimal test functions with respect [ tothequasi-optimalnormbystudyingtheirconvergencenumerically. Inordertofacilitate 1 understanding of the method, we also include a detailed explanation of the methodology v 7 from the algorithmic point of view. 4 8 Keywords: convection-diffusion; finite element method; discontinuous Petrov-Galerkin; 1 . numerical stability; robustness; automatic stabilization technique 1 0 2 1 1 Introduction : v Xi The purpose of this paper is to study the application of the discontinuous Petrov-Galerkin finite r element methodology to convection-diffusion problems of the form a (cid:15)∆u+a ∇u = f in Ω − · (1) u = g on ∂Ω whereΩisaboundedtwo-dimensionaldomain, aistheflowvelocity, (cid:15)isthediffusionparameter and f represents a source term. The function g determines the values of the solution u on the boundary of the domain denoted by ∂Ω. This equation is a basic model for various transport processes in nature and engineering. Typical applications include the movement of substances (dissolved nutrients, toxins) in the environment. Problem (1) may also be viewed as a stepping ∗[email protected] †[email protected] ‡[email protected] 1 stone for designing numerical methods for more complex problems such as the Navier-Stokes equations for compressible and incompressible flows. Construction of numerical solution schemes for (1) that can approximate convection dom- inated flows is known to be an inherently difficult task. Many scientists and engineers have addressed the problem and designed algorithms based on finite difference, finite volume, finite element and spectral approaches. In the context of the finite element method, the standard Ritz-Bubnov-Galerkin formulation is known to be deficient and various alternative formulations have been proposed in the literature, see for instance [7–10, 12, 21, 22, 27–29, 32, 34, 35, 44]. The success of the traditional Galerkin finite element method in most structural problems is based on its best approximation property. This means that the difference between the finite element solution and the exact solution is minimized with respect to a norm induced by the strain energy product. Numerical issues arise when the best approximation property is lost. This happens, for instance, when the standard Galerkin method is applied to convective transport problems. In these problems, the system matrix associated to convection is not symmetric and numerical solutions tend to show spurious, non-physical oscillations unless the finite element mesh is heavily refined. It should be noted that the accuracy of the finite element method may also degrade in problems with symmetric system matrices. Thin-body problems of solidmechanicsandwavepropagationproblemsconstitutetypicalexamples. Intheseproblems, the standard Galerkin method leads to non-physical solutions suffering from the phenomena of numerical locking (small thickness) and numerical dispersion (high wave number) unless an over-refined mesh is used. To save computer resources, alternative formulations which produce reasonable approxima- tions regardless of the mesh size have been developed. In the context of convective transport problems,awell-knownmethodistheStreamlineUpwindPetrov-Galerkin(SUPG)formulation, see [9, 36]. The method belongs to a larger class of stabilized finite element methods reviewed, for example, in [12, 25]. The stabilized finite element methods are geared mainly towards the classical h-version of the finite element method although some higher order formulations have been studied in [4, 11, 24, 31]. A new, general finite element framework has been developed by Demkowicz and Gopalakr- ishnan in [15–17]. This variational framework is of discontinuous Petrov-Galerkin (DPG) type and can be utilized, for example, to symmetrize non-symmetric problems by using the concept of optimal test functions similar to the early methodology developed by Barrett and Morton in [3], see also [21, 22, 26, 38, 39]. More precisely, if a fully discontinuous finite element space is used for the test functions, as in the DPG method of Bottasso et al. [5], the optimal test functions can be approximated locally in an enriched finite element space. The same procedure can also be used to compute local a posteriori error estimates to guide automatic adaptive mesh refinement as demonstrated in [19] in the context of convection-dominated diffusion problems and in [41] in the context of shell boundary layers. Previous work in the DPG framework applied to Problem (1), has utilized adaptive mesh refinements [14, 19] to fully resolve fine scales (boundary and interior layers) in the solution. However, in engineering applications this is not always possible. In this work, we focus on the coarse-mesh accuracy of the DPG method. That is, we investigate how the effect of fine scales to the solution can be taken into account without resolving them by the trial space mesh. It should be noted that coarse-mesh accuracy is valuable even in the context of adaptive mesh refinement to avoid unnecessary mesh refinement. We complement the earlier works [16, 19] by resolving the stable test function space with respect to the optimal (quasi-optimal) test space norm. The corresponding DPG method is expected to be uniformly stable with respect to vanishing diffusion, but the local problems 2 for the test functions become singularly perturbed and therefore more difficult to solve. The main contribution of our study is to show how this problem can be approached by utilizing a carefully designed element sub-grid discretization. The present work is a continuation of our recent work [42, 43] and extends the earlier analysis in one space dimension to two dimensions. The paper is structured to be self-contained and focused on bridging the theoretical frame- work to a practical understanding. Section 2 lays down the abstract functional analytic setup of the method. Then in Section 3, we specialize the theoretical framework to Problem (1) and introduce different test space norms. In Section 4, we re-explain the DPG framework from a purely algorithmic point of view in order to highlight the practical implications of the different choices of the test space norm. We provide details on how to construct the linear algebraic sys- tems constituting the local optimal test functions. Numerical results are presented in Section 5 and the paper ends with conclusions and remarks in Section 6. 2 Automatically Stable DPG Framework 2.1 Background The starting point of the formulation is an abstract variational statement of the form: Find u such that ∈ U (w,u) = (w), w (2) B L ∀ ∈ W where the test function space and the trial function space are assumed to be Hilbert W U spaces under inner products ( , ) and ( , ) with the corresponding norms and . W U W U · · · · ||·|| ||·|| We shall assume that the spaces are real because our application, Problem (1), deals only with real-valued functions. If the spaces and are the same and induces an inner product, then the standard W U B Galerkin method delivers the best approximation of the exact solution in the corresponding norm known as the energy norm. In the case of Problem (1) the conventional weak formulation is associated with the bilinear and linear forms (cid:90) (cid:90) (v,u) = ((cid:15)∇v ∇u+v(a ∇u)) dΩ, (v) = vf dΩ (3) B · · L Ω Ω Here does not induce an inner product because the second term is not symmetric. Conse- B quently the best approximation property in the energy norm is lost together with numerical stability, see [9, 33]. This has motivated the use of the Petrov-Galerkin method, where the combination of specially engineered test functions with added stabilization terms in the varia- tionalformulation, hasimprovedthenumericalapproximationsconsiderably. TheDPGmethod discussed here restores the best approximation property by the computation of a special set of test functions for given trial functions that produce a symmetric positive-definite algebraic system when substituted in the bilinear form. The method can be explored from many dif- ferent perspectives. It may be viewed as an automatic apparatus designed for guaranteeing the famous inf-sup condition of the Babuˇska-Brezzi theory of discrete stability, see [1, 6], or it may be interpreted as a generalized least squares method, see [17]. We treat the methodology as a numerical construction which generates a true inner-product structure for any well-posed variational problem. 3 2.2 Petrov-Galerkin Method with Optimal Test Functions ThegoalistofindalinearoperatorTwhichmapsthetrialfunctionspace tothetestfunction U space such that W (Tv,u) = (v,u) v,u (4) B A ∀ ∈ U where (v,u) is a symmetric coercive bilinear form on . If e ,e ,...,e denotes a set 1 2 n A U ×U { } of trial functions that spans a subspace of , then the use of the test functions Te in the n i U U Petrov-Galerkin method (Te ,u ) = (Te ), i = 1,...,n (5) i n i B L defines u as the orthogonal projection of the exact solution u onto the trial space with n n U respect to the inner product ( , ): A · · (v ,u u ) = 0 v n n n n A − ∀ ∈ U In other words, u is the best approximation of u in in the corresponding ‘energy’ norm, n n U that is, u u u v v n U U n ||| − ||| ≤ ||| − ||| ∀ ∈ U where (cid:112) u = (u,u) u U ||| ||| A ∀ ∈ U A canonical way to define the trial-to-test operator T : is to define Tu through U → W ∈ W the variational equality (w,Tu) = (w,u) w (6) W B ∀ ∈ W so that we have (Tv,u) = (Tv,Tu) W B and thus achieve our goal with (v,u) = (Tv,Tu) , u = Tu (7) W U W A ||| ||| || || 2.3 Localization Foreachtrialfunctione ,i = 1,...,n,Equation(6)constitutesanauxiliary,infinite-dimensional i variational problem associated with the test function space . The standard Galerkin method W maybeusedtoapproximatesolutionstoEquation(6). However, fortypicalvariationalformula- tions, like for the one corresponding to the forms of (3), Equation (6) is from the computational perspective at least as ‘heavy’ as the original problem. ThenewinsightofDemkowiczandGopalakrishnan,see[16],wastoformulatethevariational problem (2) in an ultra-weak hybrid form in such a way that the test function space becomes W fully discontinuous on a partitioning of Ω into elements K ,K ,...,K : 1 2 N { } = w : w (K) for every element K K W { | ∈ W } Aslongastheinnerproduct( , ) doesnotcoupletestfunctionsdefinedonindividualelements, W · · the auxiliary problems defined by Equation (6) can be solved locally at the element level. The localized version of the problem can be written in terms of the broken forms (cid:88) (cid:88) ( , ) = ( , )2 , ( , ) = ( , ) · · W · · W(K) B · · BK · · K K 4 which allows the global problem (6) to be written as finding Tu (K) such that ∈ W (w,Tu) = (w,u) w (K) (8) W(K) K B ∀ ∈ W This problem may then be solved approximately using the standard Galerkin method in an ˜ enriched finite element space (K) (K), that is, find Tu (K) such that n˜ n˜ W ⊂ W ∈ W ˜ (w,Tu) = (w,u) w (K) (9) W(K) K n˜ B ∀ ∈ W 2.4 Well-posedness and Robustness The above construction can be carried out provided that the original variational problem (2) is well-posed, that is, provided that the bilinear form (w,u) satisfies the conditions of the B Babuˇska-Necas-Nirenberg Theorem, which are (see [1] and [6, 13]) (w,u) C w u w , u (10) 1 W U B ≤ || || || || ∀ ∈ W ∈ U (w,u) sup B C w (11) 2 W u ≥ || || u∈U U || || (w,u) sup B C u (12) 3 U w ≥ || || w∈W W || || ThefirstconditionimpliesthattherighthandsideofEquation(6)isaboundedlinearfunctional in for every u . Therefore Tu exists and is unique. W ∈ U ∈ W The second condition guarantees that the mapping T is surjective, that is, every w ∈ W can be expressed in the form w = Tu for some u . Namely, if this was not the case, there ∈ U would exist a non-zero w˜ such that (w˜,Tu) = 0 for all u . But then also W ∈ W ∈ U (w˜,u) = 0 u B ∀ ∈ U by Equation (6) so that Equation (11) implies w˜ = 0, which is a contradiction. Finally, the coercivity of the bilinear form follows from condition (12). To see this, we A recall that the norm of an element in a Hilbert space can be computed by duality as H (v,u) H u = sup (13) H || || v v∈H H || || so that the energy norm can be expressed as (w,Tu) (w,u) W u = sup = sup B U ||| ||| w w w∈W W w∈W W || || || || Now Equation (12) implies that (u,u) C2 u 2 u A ≥ 3|| ||U ∀ ∈ U The characteristics of the energy norm depend on the values of the constants C and C 1 3 which in turn typically depend on the physical parameters and the problem geometry. As a result the Petrov-Galerkin method (5) may not be robust, that is, uniformly accurate with respect to the parameter values of a particular physical problem that is being modeled, see [2]. However, the parametric dependence can be removed completely from the bilinear form ( , ) A · · 5 by using a special norm for the test function space in (6). This norm is referred to as the optimal test space norm in the DPG literature, see for instance [48], and can be defined as w = T∗w (14) W U ||| ||| || || where T∗ : is a linear mapping such that W → U (T∗w,u) = (w,u) u U B ∀ ∈ U Equation (14) defines a proper norm in again under the assumption that the variational W problem is well-posed. Namely, the utilization of the duality argument (13), and conditions (10)–(11) show that C w T∗w C w 2 W U 1 W || || ≤ || || ≤ || || that is, is an equivalent norm to . W W |||·||| ||·|| Finally, condition (12) implies that T∗ is surjective so that a subsequent application of the duality argument (13) yields (w,u) (T∗w,u) (v,u) U U u = sup B = sup = sup = u (15) ||| |||U w T∗w v || ||U w∈W W w∈W U v∈U U ||| ||| || || || || Thus, the introduction of (14) allows us to obtain a problem for each trial function which yields a robust method (convergent in the norm as stated in (15)). This guaranteed best U ||·|| approximation property is what leads us to call the norm induced by Equation (14) optimal. As we will observe in the proceeding Section, the ultra-weak approach allows the expression of the optimal test space norm to be deduced directly from the bilinear form. Additionally, the leading terms can be expressed in a closed, localizable form. 3 DPG Method for the Convection-Diffusion Problem In this section we follow [16, 18] and derive an explicit DPG variational formulation of the convection-diffusion equation. The ultra-weak DPG variational form of the problem is associ- ated initially to a partitioning of Ω into elements denoted by Ω . We will refer to the collection h of all edges in the mesh as ∂Ω = ∂K. We will use the standard notation where L (Ω) h K 2 ∪ and L (Ω) denote the spaces of square-integrable scalar- and vector-valued functions over Ω, 2 respectively. Moreover, H1(Ω) and H(div,Ω) stand for the subspaces consisting of functions with square-integrable derivatives: H1(Ω) = v L (Ω) : ∇v L (Ω) 2 2 { ∈ ∈ } H(div,Ω) = τ L (Ω) : ∇ τ L (Ω) 2 2 { ∈ · ∈ } 3.1 Derivation of the Hybrid Ultra-Weak Variational Form The procedure begins by rewriting the second-order partial differential equation in (1) as a pair of first-order equations (cid:15)−1σ = ∇u − (16) ∇ σ +a ∇u = f · · Here we have chosen to define σ as the diffusive flux in order to be consistent with the above cited references. Other option would be to include the advective flux au in the first equation defining σ as in our earlier works [42, 43]. In our experience, the difference between the two 6 approaches appears to be mostly cosmetic. In this work we will assume that the ambient fluid is incompressible so that ∇ a = 0 and that the problem is formulated in a dimensionless form · scaled such that a 1. | | ∼ Upon multiplying the first and second equations in (16) by a vector-valued test function τ and a scalar-valued test function v, respectively, and integrating by parts element-wise, we obtain the following bilinear and linear forms corresponding to the abstract variational statement (2) (w,u) = ((cid:15)−1τ ∇v,σ) (∇ τ +a ∇v,u) + τ n,uˆ + σˆ ,v B − Ωh − · · Ωh (cid:104) · (cid:105)∂Ωh (cid:104) n (cid:105)∂Ωh (17) (w) = (v,f) L Ωh Herendenotestheoutwardunitnormaloneachelementedgeandσˆ correspondstothe(outer) n normal component of the total flux (σ + au) n. Moreover, the test function is denoted by · w = (τ,v) and the trial function by u = (σ,u,uˆ,σˆ ). For smooth functions, the notations n ( , ) and , stand for the element-wise computed inner products · · Ωh (cid:104)· ·(cid:105)∂Ωh (cid:90) (cid:90) (cid:88) (cid:88) (f,g) = fgdΩ, f,g = fgds Ωh (cid:104) (cid:105)∂Ωh K∈Ωh K K∈Ωh ∂K The element integrals ( , ) in Equation (17) make sense provided that σ,u,τ,v as well · · Ωh as ∇ τ and ∇v are square integrable over each K. Thus, the test function space is defined as · the broken space = H(div,Ω ) H1(Ω ) (18) h h W × where H(div,Ω ) = τ L (Ω) : ∇ τ L (K) K Ω h 2 2 h { ∈ · ∈ ∀ ∈ } H1(Ω ) = v L (Ω) : ∇v L (K) K Ω h 2 2 h { ∈ ∈ ∀ ∈ } Inthegeneralcase, theinterfaceaction , in(17)mustbeinterpretedasanappropriate (cid:104)· ·(cid:105)∂Ωh duality pairing between suitable generalized function spaces defined on ∂K. The right setting can be deduced from the theory of trace operators naturally associated to the spaces H(div,K) and H1(K). These spaces are formally denoted by H−1/2(∂K) and H1/2(∂K) for the normal trace τ n and the standard trace v , respectively. ∂K ∂K · | | This leads to the definition of the trial function space as = L2(Ω) L2(Ω) H1/2(∂Ω ) H−1/2(∂Ω ) (19) U × × D h × h where the spaces for the interface variables uˆ and σˆ are defined as n H1/2(∂Ω ) = v : v H1(Ω) D h { |∂Ωh ∈ D } H−1/2(∂Ω ) = η : η H(div,Ω) h { |∂Ωh ∈ } Here H1(Ω) refers to the subspace of H1(Ω) encompassing the Dirichlet boundary condition D u = g on ∂Ω. The topology of the numerical trace spaces can be characterized for instance by using the extensions with minimal norm: uˆ = inf v : v H1(Ω) such that v = uˆ || ||H1/2(∂Ωh) { || ||H1(Ω) ∈ |∂Ωh } σˆ = inf η : η H(div,Ω) such that η n = σˆ || n||H−1/2(∂Ωh) { || ||H(div,Ω) ∈ · |∂Ωh n} More details regarding these definitions can be found from [18]. For our purposes it suffices 1/2 to note that piecewise polynomial functions in H (∂Ω ) must be globally continuous on ∂Ω D h h whereas discontinuities at the vertices of Ω are allowed in the space H−1/2(∂Ω ). h h 7 3.2 Finite Element Trial Spaces We refer to P as the set of polynomials of degree lower or equal to t and to the corresponding t tensor product family by Q = P P and set Q = Q . The trial space is defined then px,py px × py p p,p as = S2 S M N Un h × h × h × h where S = v L (Ω) : v Q for every element K h 2 K p { ∈ | ∈ } M = v H1(Ω) : v P for every edge E h { ∈ D |E ∈ p+1 } N = v : v P for every edge E h E p { | ∈ } In this work we use the Bernstein polynomials as a basis for P . The t+1 basis polynomials t are defined as (cid:18) (cid:19) t B (x) = xi(1 x)t−i, i = 0,1,...,t i,t i − 3.3 Localization in the Convection-Diffusion Problem Upon selection of a norm 2 = ( , ) for the broken test function space defined in (18), a ||·||W · · W W precise form of the localized auxiliary problem (8) is obtained. This problem is solved approxi- matelyusingthestandardGalerkinmethodandanenrichedfiniteelementspacedefinedoneach ˜ ˜ element K. We will associate this space to a sub-partitioning of K into elements K ,...,K . { 1 N˜} More precisely, we define an H(div,K) H1(K)-conforming finite element space for the × test functions as (K) = Wn˜ Th˜ ×Vh˜ where = τ H(div,K) : τ Q Q for every sub-element K˜ Th˜ { ∈ |K˜∈ p˜+1,p˜× p˜,p˜+1 } = v H1(K) : v Q for every sub-element K˜ Vh˜ { ∈ |K˜∈ p˜ } Here we denote p˜= p+∆p so that the level of enrichment is controlled by both the number of ˜ sub-elements N and the increment in the order of approximation ∆p. In what follows, we will describe three different choices of the test space norm and discuss appropriate enrichments for resolving the corresponding test functions. The Standard Test Space Norm A norm known as the standard norm, or mathematician’s norm, may be constructed by com- posing the natural norms of the spaces H(div,Ω ) and H1(Ω ) as in h h w 2 = ∇ τ 2 + τ 2 + ∇v 2 + v 2 || ||SN || · ||Ωh || ||Ωh || ||Ωh || ||Ωh Thischoiceleadstotwodecoupled, standardvariationalproblemsforτ andv overeachelement K which can be solved sufficiently accurately by using an enriched finite element space with ˜ ∆p = 2 and N = 1, see [16]. The Weighted Test Space Norm While the standard norm provides a local problem that is easy to solve, it leads to approxi- mations of u that underestimate the true values considerably. It was determined in [19], that 8 the use of a special weighting of the norm near the inflow boundary could improve the accu- racy. The original version of the weighted norm, as introduced in [19], is based directly on the standard norm and reads w 2 = ∇ τ 2 + τ 2 + ∇v 2 + v 2 || ||WN || · ||Ωh,β || ||Ωh,β || ||Ωh,β || ||Ωh,β where β is a piecewise constant norm weight such that (cid:40) γ, if d(x,Γ−) δ & d(x,Γ+) δ β(x) = ≤ ≥ 1, otherwise where d(x,Γ−) and d(x,Γ+) stand for the (normal) distances between the point x and the inflow and outflow/no-flow boundaries, respectively. These are defined as Γ− = x ∂Ω : a n < 0 { ∈ · } Γ+ = x ∂Ω : a n 0 { ∈ · ≥ } It should be noted that different refined variants of the weighted norm have been recently introduced and studied in [20]. These refinements have resemblance to the quasi-optimal norm introduced next. The Quasi-Optimal Test Space Norm While the weighted test space norm resolves some of the deficiencies of the standard norm, the approach is not optimal in the sense of the DPG theory as outlined in Section 2.4. In contrast to the weighted test space norm, we describe here the optimal test space norm (14). However, the true optimal norm leads to a series of global problems so that some numerical modifications are necessary to make the approach suitable for practical computations. This way we arrive to a quasi-optimal test space norm which retains the essential features of the optimal test space norm but is localizable. Let the trial space norm be the natural norm arising from the definition (19) (cid:110) (cid:111)1/2 u = σ 2 + u 2 + uˆ 2 + σˆ 2 || ||U || ||L2(Ω) || ||L2(Ω) || ||H1/2(∂Ωh) || n||H−1/2(∂Ωh) Thus, using the definition (17) of the bilinear form, we can express the optimal test space norm (14) as w 2 = (cid:15)−1τ ∇v 2 + ∇ τ +a ∇v 2 + [τ n] 2 + [vn] 2 ||| |||W || − ||Ωh || · · ||Ωh || · ||∂Ωh || ||∂Ωh where uˆ,τ n v,σˆ [τ n] = sup (cid:104) · (cid:105)∂Ωh, [vn] = sup (cid:104) n(cid:105)∂Ωh || · ||∂Ωh uˆ || ||∂Ωh σˆ uˆ∈HD1/2(∂Ωh) || ||H1/2(∂Ωh) σˆn∈H−1/2(∂Ωh) || n||H−1/2(∂Ωh) The last two terms can be made more precise, see [18], but the important point is that they involve the ‘jumps’ of τ and v along the inter-element boundaries and therefore prevent the element-wise computation of the optimal test functions. However, a localizable variant of the optimalnormmayobtainedsimplybydroppingthesetermsandaddinganintegralcontribution of v, so as to remove the ‘zero-energy mode’ (τ,v) = (0,const.). This leads us to the so-called quasi-optimal test space norm w 2 = (cid:15)−1τ ∇v 2 + ∇ τ +a ∇v 2 +α τ 2 +α v 2 (20) || ||QON || − ||Ωh || · · ||Ωh 1|| ||Ωh 2|| ||Ωh 9 where α 0 and α > 0 are numerical regularization parameters to be chosen. Notice that 1 2 ≥ while the addition of a contribution of τ 2 is not necessary, we have observed that a proper || ||Ωh regularization of τ improves the accuracy of the final approximation. In our algorithm we choose α = (cid:15)−3/2 and α = 1. 1 2 A critical issue in the application of the quasi-optimal norm (20) is that the auxiliary problems for the test functions become singularly perturbed. This was explored in the one- dimensionalcasein[42,43]whereaspeciallydesignedsub-gridmeshwasemployedtoresolvethe potential boundary layers in the test functions. To extend the idea to two spatial dimensions, we observe that the strong form of the problem (8) associated to the quasi-optimal test space norm (20) takes the form (cid:40) ∇(∇ τ)+((cid:15)−2 +α )τ ∇(a ∇v) (cid:15)−1∇v = ε−1σ +∇u 1 − · − · − (21) ∆v a ∇(a ∇v)+α v a ∇(∇ τ)+(cid:15)−1∇ τ = ∇ σ +a ∇u 2 − − · · − · · · · · in every K. Notice that the system (21) corresponds to the Euler-Lagrange equations for the variational problem (6) and is accompanied with inhomogeneous natural boundary conditions involving also the interface variables σˆ and uˆ. n While a detailed regularity and asymptotic analysis of the system (21) is out of the scope of this paper, a preliminary examination of the system when K = (0,1) (0,1) and a = (1,0) or × a = (0,1) reveals that the solution may exhibit regular, exponential boundary layers behaving like e−λx1, e−λ(1−x1), e−λx2, e−λ(1−x2) where the characteristic exponent is given by λ = (cid:15)−1√1+α (cid:15)2. Notice that when α = (cid:15)3/2, 1 1 we have λ = (cid:15)−1+ ((cid:15)−1/2). To account for these features in the automatic computation of the O optimal test functions, we employ a layer-adapted Shishkin mesh on each K as shown in Figure 1. Such meshes are commonly used in the approximation of singularly perturbed problems, see [37, 46]. We follow a strategy which consists of adding one layer of ’needle’ elements of width (p˜(cid:15)) near the boundary and has been proven to be optimal in the context of the hp-version O of FEM for reaction-diffusion equations, see [40, 45, 47]. We emphasize that the sub-grid discretization is needed only when using the quasi-optimal test space norm. Failure to address the fine-scale features of the corresponding test functions (here using the sub-grid) degrades the accuracy of the final solution and dissipates all benefits of the quasi-optimal test norm. 4 Algorithmic Considerations At the high level, the DPG method with optimal test functions is similar to the standard finite element method detailed in Algorithm 1. Initially, the mesh information is input and mappings from the local stiffness matrix to the global matrix are determined. After this, we loop over the elements and form the element contributions to the global stiffness matrix and load vector. The main conceptual difference between the DPG method and standard FEM is the following: the local stiffness matrix will be computed as a by-product of computing the optimal test functions. Another difference is that, because of the hybridization, the trial spaces are non-standard and require some specialized handling. The remaining portion of this section will describe both differences in more detail. For a review of standard FEM implementation, see [14, 30]. 10