Stability and Functional Superconvergence of Narrow-Stencil Second-Derivative Generalized Summation-By-Parts Discretizations
SStability and Functional Superconvergence of Narrow-StencilSecond-Derivative Generalized Summation-By-Parts Discretizations
Zelalem Arega Worku * David W. Zingg * Abstract
We analyze the stability and functional superconvergence of discretizations of diffusion problems withthe narrow-stencil second-derivative generalized summation-by-parts (SBP) operators coupled with simultaneousapproximation terms (SATs). Provided that the primal and adjoint solutions are sufficiently smooth and theSBP-SAT discretization is primal and adjoint consistent, we show that linear functionals associated with thesteady diffusion problem superconverge at a rate of 2 p when a degree p + 1 narrow-stencil or a degree p wide-stencil generalized SBP operator is used for the spatial discretization. Sufficient conditions for stability of adjointconsistent discretizations with the narrow-stencil generalized SBP operators are presented. The stability analysisassumes nullspace consistency of the second-derivative operator and the invertibility of the matrix approximatingthe first derivative at the element boundaries. The theoretical results are verified by numerical experiments withthe one-dimensional Poisson problem. Keywords
Summation-by-parts · Adjoint consistency · Simultaneous approximation term · Narrow-stencil · Functional superconvergence
Mathematics Subject Classification (2010) · · · Compared to wide-stencil summation-by-parts (SBP) operators, explicitly formed narrow-stencil second-derivativeSBP operators provide smaller solution error, superior solution convergence rates, compact stencil width, and bet-ter damping of high frequency modes [29,30,26,12,13]. As with the wide-stencil operators, narrow-stencil second-derivative operators are coupled by simultaneous approximation terms (SATs) [6]. However, the SAT coefficientsderived for wide-stencil SBP operators must be modified for implementations with narrow-stencil SBP operatorsto achieve stability and adjoint consistency simultaneously. Unfortunately, the analysis required to find such SATcoefficients for narrow-stencil SBP operators is more involved, e.g., see [13].Hicken and Zingg [22] showed that adjoint consistent SBP-SAT discretizations of linear elliptic partial differ-ential equations (PDEs) lead to functional superconvergence (see also [5,24,20]). In their study, they analyzeddiscretizations with wide-stencil second-derivative classical SBP (CSBP) operators by posing second-order linearPDEs as a system of first-order equations and determined the conditions that the SATs must satisfy for ad-joint consistency and functional superconvergence. A similar analysis is conducted in [36] for multidimensionalSBP operators, but without posing the second-order linear PDEs as a system of first-order equations. The latterapproach enables analysis of functional accuracy of adjoint consistent discretizations of diffusion problems withnarrow-stencil SBP operators. While the stability of discretizations arising from narrow-stencil second-derivativeoperators is well-studied (e.g., see [7,30,26,16,28,27]), it is only recently (see, e.g., [13,14]) that conditions for which such discretizations satisfy both stability and adjoint consistency requirements are presented. Eriksson [13],used the eigendecomposition technique to find the conditions on the SAT coefficients that enable constructionof stable and adjoint consistent discretizations of diffusion problems. In a subsequent paper [14], Eriksson andNordstr¨om used a variant of the approach in [13] to find a more general set of SAT coefficients. Although theseSAT coefficients lead to adjoint consistent and stable discretizations in practice, the analysis in [14] assumes acondition that is not satisfied by many narrow-stencil second-derivative operators in the literature, including thosein [29,28,26,12,27]. Furthermore, it is not straightforward how the theory extends to narrow-stencil generalizedSBP operators which have one or more of the following characteristics: exclusion of one or both boundary nodes,non-repeating interior point operators, and non-uniform nodal distribution [12]. Zelalem Arega WorkuE-mail: [email protected] W. ZinggE-mail: [email protected] * Institute for Aerospace Studies, University of Toronto, Toronto, Ontario, M3H 5T6, Canada Second-derivative operators formed by applying first-derivative operators twice. Also known as compact-stencil second-derivative operators. a r X i v : . [ m a t h . NA ] F e b Z. Worku, D.W. Zingg
The first objective of this paper is to establish the conditions required for the stability of adjoint consistent SBP-SAT discretizations of diffusion problems with the generalized narrow-stencil second-derivative SBP operators. Weuse the “borrowing trick”[7] in the energy stability analysis which directly applies to the diagonal- and block-norm narrow-stencil second-derivative SBP operators in [29,28,26,27] and to the generalized SBP operators of Del ReyFern´andez and Zingg [12] upon minor modifications of the derivative operators at element boundaries. The secondobjective is to show that primal and adjoint consistent discretizations lead to functional convergence rates of 2 p when a degree p + 1 narrow-stencil or a degree p wide-stencil diagonal-norm second-derivative generalized SBPoperator is used to discretize steady diffusion problems for which the primal and adjoint solutions are sufficientlysmooth. We also show that the functional converges at a rate of 2 p irrespective of whether or not the schemeis adjoint consistent when a degree 2 p − We closely follow the notation used in [12,10,37,36]. A one-dimensional compact domain is considered, and it istessellated into n e non-overlapping elements, T h := {{ Ω k } n e k =1 : Ω = ∪ n e k =1 Ω k } . The boundaries of each element willbe referred to as interfaces, and we denote their union by Γ k := ∂Ω k . The set of all interior interfaces is denoted by Γ I := { Γ k ∩ Γ v : k, v = 1 , . . . , n e , k (cid:54) = v } , while the element interfaces for which Dirichlet and Neumann boundaryconditions are enforced are in the sets Γ D and Γ N , respectively, and Γ := Γ I ∪ Γ D ∪ Γ N . Operators associated withthe left and right interfaces of Ω k bear the subscripts (cid:96) and r , respectively, and the left and right most elementsare indicated by the subscripts L and R , respectively, e.g., D (cid:96)L is a derivative operator at the left interface of theleft most element. The set of n p volume nodes in element Ω k is represented by x k = { x i } n p i =1 . Uppercase scripttype, e.g., U k ∈ C ∞ ( Ω k ), is used for continuous functions, and P p ( ˆ Ω ) denotes the space of polynomials up to totaldegree p , which has a cardinality of n ∗ p = p + 1. Bold letters, e.g., u k ∈ R n p , delineate the restriction of U k to gridpoints x k , while solution vectors to the discrete systems of equations have subscript h , e.g., u h,k ∈ R n p . For thepurpose of the functional convergence analysis in Section 4.3, we define h := max a,b ∈ x k | a − b | as the size of anelement. Matrices are denoted by sans-serif uppercase letters, e.g., V ∈ R n p × n p ; denotes a vector consisting ofall ones, denotes a vector or matrix consisting of all zeros. The sizes of and should be clear from context.Definitions of the first- and second-derivative SBP operators presented in [12] are stated below. For the con-struction of narrow-stencil second-derivative SBP operators, we refer the reader to [10,12,26,27]. Definition 2.1 (Generalized first-derivative SBP operator)
The matrix D k ∈ R n p × n p is a degree p SBPoperator approximating the first derivative ∂∂x on the set of nodes x k , which need neither be uniform nor includenodes on the boundaries and may have nodes outside the domain of element Ω k , if [12]1. D k p = ∂ P ∂x for all P ∈ P p ( Ω k )2. D k = H − k Q k , where H k is a symmetric positive definite (SPD) matrix, and3. Q k = S k + E k , where S k = − S Tk , E k = E Tk , and E k satisfies p T E k q = (cid:80) γ ∈ Γ k [ P ] γ [ Q ] γ n γk for all P , Q ∈ P τ ( Ω k ),where τ ≥ p , and n γk = 1 if γ is the right interface of Ω k , otherwise n γk = − H k , may be diagonal or dense. A dense-norm matrix refers to any norm matrix that is notdiagonal, which includes the block-norm matrix. The block-norm matrix has diagonal entries at the interior points(containing h ) and dense blocks at the top-left and bottom-right corners corresponding to the boundary nodes.The L inner product of two functions P and Q is approximated by [23,10,21,11] p T H k q = (cid:90) Ω k PQ d Ω + O (cid:16) h p (cid:17) , and H k defines the norm u T H k u = (cid:107) u (cid:107) H = (cid:90) Ω k U d Ω + O (cid:16) h p (cid:17) . Also referred to as full-norm matrix.tability and Functional Superconvergence of Narrow-Stencil SBP Discretizations 3
The E k matrix is constructed as [10,11] E k = (cid:88) γ ⊂ Γ k n γk R Tγk R γk = R Trk R rk − R T(cid:96)k R (cid:96)k , (2.1)where R γk is an extrapolation row vector of at least order h τ +1 accuracy, i.e., R γk u k = [ U k ] γ + O ( h ≥ τ +1 ).Furthermore, we define an operator that extrapolates the product of the diffusion coefficient and the derivative ofthe solution from volume nodes to an interface as D γk = n γk R γk Λ k D b,k . (2.2) Definition 2.2 (Order-matched narrow-stencil second-derivative generalized SBP operator)
Thenarrow-stencil second-derivative operator D (2) k of degree p + 1, approximating ∂∂x ( λ k ∂ U k ∂x ), is order-matched withthe first-derivative operator D k = H − k Q k of degree p on the nodal set x k if [12] D (2) k ( λ k ) p k = ∂∂x (cid:18) λ k ∂ P k ∂x (cid:19) , ∀ ( λ k P k ) ∈ P p +1 ( Ω k ) , (2.3)and D (2) k is of the form D (2) k = H − k [ − M k + E k Λ k D b,k ] , (2.4)where M k = (cid:80) n p i =1 Λ k ( i, i ) ¯ M i , ¯ M i are symmetric positive semidefinite matrices, Λ k = diag( λ k ( x ) , λ k ( x ) , ..., λ k ( x n p )) , and D b,k is an approximation to the first derivative of degree and order ≥ p + 1.The order-matched SBP operators in Definition 2.2 are assumed to have a diagonal-norm matrix. Note thatfor the m th derivative, the degree and order are related by order = degree − m + 1; consequently, both thediagonal-norm narrow-stencil D (2) k and D k operators are order p accurate, while a diagonal-norm wide-stencilsecond-derivative operator, which has the decomposition D k Λ k D k = H − k [ − D Tk H k Λ k D k + E k Λ k D k ] , (2.5)is order p − p centered-difference interior operator. At the boundaries, however, the block-norm wide- and narrow-stencilsecond-derivative operators are closed with order 2 p − p − p one-sidedstencils used with the diagonal-norm wide- and narrow-stencil SBP operators, respectively. Furthermore, the D b,k matrix of a block-norm operator contains order 2 p − Remark 2.1
In this work, we do not assume that M k is necessarily symmetric positive semidefinite; rather weassume that M k + M Tk is symmetric positive semidefinite, which allows the analysis to be extended to a moregeneral class of explicitly formed second-derivative operators, including the block-norm SBP operators in [27,29], which do not have symmetric M k matrix.Another decomposition of second-derivative SBP operators, which is instrumental for the adjoint consistencyand functional superconvergence analyses in Section 4, is presented below. Proposition 2.1
A second-derivative operator of the form (2.4) , for which M k is not necessarily symmetric, canbe decomposed as D (2) k = H − k (cid:16) D (2) k (cid:17) T H k − H − k D Trk R rk − H − k D T(cid:96)k R (cid:96)k + H − k R Trk D rk + H − k R T(cid:96)k D (cid:96)k − H − k (cid:16) M k − M Tk (cid:17) . (2.6) Proof
Substituting (2.1) and (2.2) into (2.4), we have D (2) k = H − k [ − M k + E k Λ k D b,k ] = − H − k M k + H − k (cid:16) R Trk R rk − R T(cid:96)k R (cid:96)k (cid:17) Λ k D b,k = − H − k M k + H − k R Trk D rk + H − k R T(cid:96)k D (cid:96)k Z. Worku, D.W. Zingg
Adding and subtracting H − k (cid:16) D (2) k (cid:17) T H k , we get D (2) k = − H − k M k + H − k (cid:16) D (2) k (cid:17) T H k − H − k (cid:16) D (2) k (cid:17) T H k + H − k R Trk D rk + H − k R T(cid:96)k D (cid:96)k = H − k (cid:16) D (2) k (cid:17) T H k − H − k (cid:104) − M Tk + D Tb,k Λ Tk E Tk (cid:105) + H − k R Trk D rk + H − k R T(cid:96)k D (cid:96)k − H − k M k = H − k (cid:16) D (2) k (cid:17) T H k − H − k D Trk R rk − H − k D T(cid:96)k R (cid:96)k + H − k R Trk D rk + H − k R T(cid:96)k D (cid:96)k − H − k (cid:16) M k − M Tk (cid:17) , which is the desired result. (cid:117)(cid:116) We consider the one-dimensional diffusion problem ∂ U ∂t − ∂∂x (cid:18) λ ∂ U ∂x (cid:19) = F ∀ x ∈ Ω, U = U at t = 0 , U| Γ D = U D , n γ (cid:18) λ ∂ U ∂x (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) Γ N = U N , (3.1)where F ∈ L ( Ω ), λ = λ ( x ) is a positive diffusivity coefficient, and Γ D is not empty. For functional error analysisand numerical experiment purposes, we consider the steady version of (3.1), the Poisson problem. We also considera compatible linear functional of the form I ( U ) = (cid:90) Ω GU d Ω − ψ D (cid:20) λ ∂ U ∂x n γ (cid:21) Γ D + ψ N U| Γ N , (3.2)where G ∈ L ( Ω ), ψ N = n γ ( λ ∂ψ∂x ) ∈ L ( Γ N ), and ψ D ∈ L ( Γ D ). A linear functional is compatible with the steadyversion of (3.1) if [18] (cid:90) Ω ψ ∂∂x (cid:18) λ ∂ U ∂x (cid:19) d Ω + U D (cid:20) λ ∂ψ∂x n γ (cid:21) Γ D − U N ψ | Γ N = (cid:90) Ω U ∂∂x (cid:18) λ ∂ψ∂x (cid:19) d Ω + ψ D (cid:20) λ ∂ U ∂x n γ (cid:21) Γ D − ψ N U| Γ N , (3.3)i.e., I ( U ) = I ( ψ ) = (cid:90) Ω ψ F d Ω − U D (cid:20) λ ∂ψ∂x n γ (cid:21) Γ D + U N ψ | Γ N . (3.4)Under the compatibility condition on the functional, the adjoint, ψ , satisfies the PDE (see, e.g., [22,37,18]) − ∂∂x (cid:18) λ ∂ψ∂x (cid:19) = G ∀ x ∈ Ω, ψ | Γ D = ψ D , n γ (cid:18) λ ∂ψ∂x (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) Γ N = ψ N . (3.5)The SBP-SAT semi-discretization of the diffusion problem, (3.1), is given byd u h,k d t = D (2) k u h,k + f k − H − k s Ik ( u h,k ) − H − k s Bk ( u h,k , u D , u N ) =: R h,u , (3.6)where f k is the restriction of F to the volume nodes in Ω k and the interface SATs, s Ik , and boundary SATs, s Bk , given in [37,36] are specialized for one-dimensional implementation as s Ik ( u h,k ) = (cid:88) γ ⊂ Γ Ik (cid:2) R Tγk D Tγk (cid:3) (cid:34) T (1) γk T (3) γk T (2) γk T (4) γk (cid:35) (cid:20) R γk u h,k − R γv u h,v D γk u h,k + D γv u h,v (cid:21) (3.7)and s Bk ( u h,k , u D , u N ) = (cid:40)(cid:2) R Tγk D Tγk (cid:3) (cid:34) T ( D ) γk − (cid:35) ( R γk u h,k − u D ) (cid:41) γ ⊂ Γ D + (cid:110) R Tγk ( D γk u h,k − u N ) (cid:111) γ ⊂ Γ N . (3.8)The SAT coefficients T (1) γk , T (2) γk , T (3) γk , T (4) γk , T ( D ) γk ∈ R are determined such that the scheme satisfies desired propertiessuch as conservation, adjoint consistency, and energy stability. For implementations with wide-stencil operators,we replace D (2) k by D k Λ k D k in (3.6) and D b,k by D k in (2.2).Substituting the restriction of the exact solution to grid points, u k , in (3.6) to (3.8), we see that the right-handside (RHS) of (3.6) yields a discretization error of O ( h p ) when an order-matched narrow-stencil second-derivativeSBP operator is used; hence, the discretization of the primal problem is consistent. In contrast, for diagonal-norm wide-stencil SBP operators, the discretization error is O ( h p − ) while for block-norm second-derivative SBP operators, it reduces to O ( h p − ). tability and Functional Superconvergence of Narrow-Stencil SBP Discretizations 5 In this section, we present the two main results of this paper. After establishing the conditions required for adjointconsistency and conservation, we show that primal and adjoint consistent SBP-SAT discretizations of the Poissonproblem with the diagonal-norm narrow-stencil second-derivative operators lead to functional superconvergence.To achieve this goal, we closely follow the technique used to show functional superconvergence in [36]. Then, weuse the energy method to find sufficient conditions that the SATs must satisfy for the stability of discretizationswith narrow-stencil generalized SBP operators before stating a few concrete examples of such SATs.4.1 Adjoint ConsistencyAdjoint consistency requires that the discrete adjoint problem, (cid:88) Ω k ∈T h (cid:0) L ∗ h,k ( ψ h ) − g k (cid:1) = , (4.1)where L ∗ h,k is the discrete adjoint operator, corresponding to the steady version of the primal problem (3.6) satisfylim h → (cid:88) Ω k ∈T h (cid:13)(cid:13) L ∗ h,k ( ψ k ) − g k (cid:13)(cid:13) H k = 0 . (4.2)To find the discrete adjoint operator, we begin by discretizing the two forms of the functional, (3.2) and (3.4),as I h ( u h ) = (cid:88) Ω k ∈T h g Tk H k u h,k − ψ D D (cid:96)L u h,L + ψ N R rR u h,R + ψ D T ( D ) (cid:96)L ( R (cid:96)L u h,L − u D ) , (4.3) I h ( ψ h ) = (cid:88) Ω k ∈T h f Tk H k ψ h,k − u D D (cid:96)L ψ h,L + u N R rR ψ h,R + u D T ( D ) (cid:96)L ( R (cid:96)L ψ h,L − ψ D ) , (4.4)where we have assumed that the Dirichlet and Neumann boundary conditions are enforced on the left and rightboundaries, respectively. The last terms in (4.3) and (4.4) arise from consistent modifications of the functional,see [18,22,20,37,36]. Note that in cases where a Dirichlet boundary condition is enforced on both boundaries, weapply the Dirichlet SATs given in (3.8) on both boundaries and modify the discrete functionals, (4.3) and (4.4),by replacing the Neumann boundary terms by Dirichlet right boundary terms similar to those given for the leftboundary. The theory developed holds for such cases without significant modification.To derive the conditions required for adjoint consistency, we set I h ( u h ) − I h ( ψ h ) = 0, which is a discreteanalogue of the relation I ( U ) − I ( ψ ) = 0. Adding (cid:80) Ω k ⊂T h ψ Th,k H k R h,k + I h ( ψ h ) − I h ( ψ h ) = 0 to the RHS of(4.3) and rearranging we find I h ( u h ) = I h ( ψ h ) + (cid:88) Ω k ⊂T h g Tk H k u h,k − ψ D D (cid:96)L u h,L + ψ N R rR u h,R − u D T ( D ) (cid:96)L ( R (cid:96)L ψ h,L − ψ D )+ ψ D T ( D ) (cid:96)L ( R (cid:96)L u h,L − u D ) + u D D (cid:96)L ψ h,L − u N R rR ψ h,R + (cid:88) Ω k ⊂T h (cid:104) ψ Th,k H k D (2) k u h,k − ψ Th,k s Ik ( u h,k ) − ψ Th,k s Bk ( u h,k , u D , u N ) (cid:105) . (4.5) Transposing (4.5), enforcing I h ( u h ) − I h ( ψ h ) = 0, applying identity (2.6), and simplifying, we obtain (cid:88) Ω k ⊂T h (cid:26) u Th,k H k (cid:16) D (2) k ψ h,k + g k (cid:17) + u Th,k (cid:16) M k − M Tk (cid:17) ψ h,k (cid:27) − u Th,L R T(cid:96)L T ( D ) (cid:96)L ( R (cid:96)L ψ h,L − ψ D ) − (cid:88) γ ⊂ Γ I R γk u h,k R γv u h,v D γk u h,k D γv u h,v T T (1) γk − T (1) γv T (2) γk + 1 − T (2) γv − T (1) γk T (1) γv − T (2) γk T (2) γv + 1 T (3) γk − T (3) γv T (4) γk T (4) γv T (3) γk T (3) γv − T (4) γk T (4) γv R γk ψ h,k R γv ψ h,v D γk ψ h,k D γv ψ h,v + u Th,L D T(cid:96)L ( R (cid:96)L ψ h,L − ψ D ) − u Th,R R TrR ( D rR ψ h,k − ψ N ) = 0 , (4.6)from which we extract the discrete adjoint operator on element Ω k as L ∗ h,k ( ψ h ) = − D (2) k ψ h,k − H − k ( M k − M Tk ) ψ h,k + H − k ( s Ik ) ∗ ( ψ h,k ) + H − k ( s Bk ) ∗ ( ψ h,k , ψ D , ψ N ) , (4.7) Z. Worku, D.W. Zingg where the interface and boundary SATs for the adjoint problem are given, respectively, by (cid:16) s Ik (cid:17) ∗ = (cid:88) γ ⊂ Γ Ik (cid:2) R Tγk D Tγk (cid:3) (cid:34) T (1) γk − T (1) γv T (2) γk + 1 − T (2) γv T (3) γk − T (3) γv T (4) γk T (4) γv (cid:35) R γk ψ h,k R γv ψ h,v D γk ψ h,k D γv ψ h,v , (4.8) (cid:16) s Bk (cid:17) ∗ = (cid:40)(cid:2) R Tγk D Tγk (cid:3) (cid:34) T ( D ) γk − (cid:35) (cid:2) R γk ψ h,k − ψ D (cid:3)(cid:41) γ ⊂ Γ D + (cid:110) R Tγk ( D γk ψ h,k − ψ N ) (cid:111) γ ⊂ Γ N . (4.9)Furthermore, we define the residual of the SBP-SAT discretization of the adjoint problem as R h,ψ := D (2) k ψ h,k + g k + H − k ( M k − M Tk ) ψ h,k − H − k ( s Ik ) ∗ ( ψ h,k ) − H − k ( s Bk ) ∗ ( ψ h,k , ψ D , ψ N ) = , (4.10)Substituting the exact adjoint solution into (4.10), we observe that R h,ψ is O ( h ≥ p − ), i.e., the discretizationof the adjoint problem is consistent, if T (1) γk = T (1) γv , T (2) γk + 1 = − T (2) γv , T (3) γk − − T (3) γv , T (4) γk = T (4) γv , M k = M Tk . (4.11)For discretizations with wide-stencil second-derivative operators, the last condition, M k = M Tk , is satisfied bydefault.4.2 ConservationFor conservation, the homogeneous diffusion problem (3.1), i.e., F = 0, should satisfy Gauss’s theorem discretely,i.e., (cid:80) Ω k ⊂T h T H k d u k / d t must depend only on the boundary terms. Premultiplying R h,u defined in (3.6) by T H k , setting f k = 0, summing over all elements, and applying the decomposition of D (2) k given in (2.4) yields (cid:88) Ω k ⊂T h T H k R h,u = − (cid:88) γ ⊂ Γ I T T (1) γk − T (1) γk T (3) γk − T (3) γk − T (1) γv T (1) γv T (3) γv T (3) γv − T (2) γk − T (2) γk T (4) γk T (4) γk − T (2) γv T (2) γv T (4) γv T (4) γv R γk u h,k R γv u h,v D γk u h,k D γv u h,v − (cid:88) Ω k ⊂T h T M k u h,k − (cid:40)(cid:20) (cid:21) T (cid:20) T Dγ − − (cid:21) (cid:20) R γk u k − u D D γk u k (cid:21)(cid:41) γ ⊂ Γ D + u N , (4.12)which reduces to a sum of boundary terms only, (cid:88) Ω k ⊂T h T H k R h,u = (cid:110) D γk u k − T Dγ ( R γk u k − u D ) (cid:111) γ ⊂ Γ D + u N , (4.13)as required for conservation of the discretization if T (1) γk = T (1) γv , T (3) γk − − T (3) γv , T M k = . (4.14)Comparing (4.14) and (4.11) and noting that M k = , we see that adjoint consistency implies conservation, asnoted in [2,19,36].4.3 Functional SuperconvergenceWithout loss of generality, we assume that the domain is tessellated using two elements, Ω L and Ω R . In thesubsequent analysis, we will use the vectors u , u h , ψ , ψ h , f , g , E ( u h ), F ( ψ h ) ∈ R n p given by u h = (cid:20) u h,L u h,R (cid:21) , ψ h = (cid:20) ψ h,L ψ h,R (cid:21) , u = (cid:20) u L u R (cid:21) , ψ = (cid:20) ψ L ψ R (cid:21) , f = (cid:20) f L f R (cid:21) , g = (cid:20) g L g R (cid:21) , (4.15) E ( u h ) = (cid:34)(cid:16) R T(cid:96)L T ( D ) (cid:96)L − D T(cid:96)L (cid:17) ( R (cid:96)L u h,L − u D ) R TrR ( D rR u h,R − u N ) (cid:35) , F ( ψ h ) = (cid:34)(cid:16) R T(cid:96)L T ( D ) (cid:96)L − D T(cid:96)L (cid:17) ( R (cid:96)L ψ h,L − ψ D ) R TrR ( D rR ψ h,R − ψ N ) (cid:35) , (4.16) tability and Functional Superconvergence of Narrow-Stencil SBP Discretizations 7 and the matrices A , B , H , D (2) ∈ R n p × n p with block entries A = (cid:20) R rL D rL (cid:21) T (cid:34) T (1) rL T (3) rL T (2) rL T (4) rL (cid:35) (cid:20) R rL D rL (cid:21) , B = (cid:20) R rL D rL (cid:21) T (cid:34) T (1) rL T (2) rL + 1 T (3) rL − T (4) rL (cid:35) (cid:20) R rL D rL (cid:21) , A = (cid:20) R rL D rL (cid:21) T (cid:34) T (1) rL T (3) rL T (2) rL T (4) rL (cid:35) (cid:20) − R (cid:96)R D (cid:96)R (cid:21) , B = (cid:20) R rL D rL (cid:21) T (cid:34) T (1) rL T (2) rL + 1 T (3) rL − T (4) rL (cid:35) (cid:20) − R (cid:96)R D (cid:96)R (cid:21) , A = (cid:20) R (cid:96)R D (cid:96)R (cid:21) T (cid:34) T (1) (cid:96)R T (3) (cid:96)R T (2) (cid:96)R T (4) (cid:96)R (cid:35) (cid:20) − R rL D rL (cid:21) , B = (cid:20) R (cid:96)R D (cid:96)R (cid:21) T (cid:34) T (1) (cid:96)R T (2) (cid:96)R + 1 T (3) (cid:96)R − T (4) (cid:96)R (cid:35) (cid:20) − R rL D rL (cid:21) , A = (cid:20) R (cid:96)R D (cid:96)R (cid:21) T (cid:34) T (1) (cid:96)R T (3) (cid:96)R T (2) (cid:96)R T (4) (cid:96)R (cid:35) (cid:20) R (cid:96)R D (cid:96)R (cid:21) , B = (cid:20) R (cid:96)R D (cid:96)R (cid:21) T (cid:34) T (1) (cid:96)R T (2) (cid:96)R + 1 T (3) (cid:96)R − T (4) (cid:96)R (cid:35) (cid:20) R (cid:96)R D (cid:96)R (cid:21) , H = (cid:20) H L H R (cid:21) , D (2) = (cid:34) D (2) L D (2) R (cid:35) , M = (cid:20) M L − M TL M R − M TR (cid:21) . (4.17)Note that for adjoint consistent schemes, it can be shown, using (4.11), that A T = B , and A T = B . (4.18)The discrete residuals corresponding to the steady version of (3.1) and the adjoint problem (3.5) can now bewritten, respectively, as R h,u ( u h ) = − D (2) u h − f + H − A u h + H − E ( u h ) = , (4.19) R h,ψ ( ψ h ) = − D (2) ψ h − g + H − B ψ h + H − F ( ψ h ) − H − M ψ h = . (4.20)Before stating the main result, we present an assumption regarding the primal and adjoint solution accuracy. Assumption 1
We assume that unique numerical solutions for the steady version of the discrete primal equation (3.6) and the discrete adjoint problem (4.10) exist, and these solutions are at least order h p +1 accurate in themaximum norm, i.e., (cid:107) u − u h (cid:107) ∞ = O ( h ≥ p +1 ) and (cid:107) ψ − ψ h (cid:107) ∞ = O ( h ≥ p +1 ) . Assumption 1 is not necessary if pointwise stability of the SBP-SAT discretization for the Poisson problem canbe demonstrated, see [17,34,35,20,22,9]. Numerical experiments with adjoint consistent discretizations show aprimal and adjoint solution convergence rate of p + 1 when a degree p diagonal-norm wide-stencil second-derivativeSBP operator is used. In contrast, a primal and adjoint solution convergence rate of p + 2 is observed when adegree p + 1 order-matched narrow-stencil second-derivative SBP operator is used with adjoint consistent SATs.The block-norm wide- and narrow-stencil second-derivative operators, on the other hand, exhibit a primal solutionconvergence rate of 2 p .We present the order of accuracy of the discrete functional approximating I ( U ) = I ( ψ ) in the followingtheorem. Theorem 4.1
Let the primal solution of the steady version of (3.1) and the adjoint solution of (3.5) be U , ψ ∈C p +2 ( Ω ) , respectively, the variable coefficient in (3.1) and (3.5) be λ ∈ C p +1 ( Ω ) , and the source terms in (3.1) and (3.5) be F , G ∈ C p ( Ω ) , respectively. If u h , ψ h ∈ R n e n p are solutions to consistent discretizations of the steadyversion of (3.1) and (3.5) , respectively, and Assumption 1 holds, then the discrete functionals (4.3) and (4.4) areorder h p accurate approximations to the compatible linear functional I ( U ) = I ( ψ ) given by (3.2) and (3.4) , i.e., I ( U ) − I h ( u h ) = O (cid:16) h p (cid:17) , (4.21) I ( ψ ) − I h ( ψ h ) = O (cid:16) h p (cid:17) . (4.22) Proof
It is sufficient to show that the result holds for a domain tessellated by two elements, Ω L and Ω R , asthe interface SATs considered couple immediate neighboring elements only. We let the Dirichlet and Neumannboundary conditions be implemented at the left and right boundaries of the domain. The boundary terms in bothforms of the functional, (3.2) and (3.4), involve the products ψλ∂ U /∂x and U λ∂ψ/∂x . Using the continuity of ψ , U , and λ , we can approximate ( ψλ∂ U /∂x ) ∈ C p +1 ( Ω ) and ( U λ∂ψ/∂x ) ∈ C p +1 ( Ω ) at the boundary nodes bydegree ≤ p polynomials. The integrands in the volume integrals of (3.2) and (3.4) are 2 p times differentiable,i.e., ( GU ) , ( ψ F ) ∈ C p ( Ω ). Since integrals are approximated by quadratures of order h p , replacing ( GU ) , ( ψ F ) ∈C p ( Ω ) by ( (cid:102) GU ) , ( (cid:103) ψ F ) ∈ P p − ( Ω ) in the functionals introduces an error of order h p . Therefore, we consider (cid:101) U , ( (cid:94) λ∂ U /∂x ) ∈ P p ( Ω ) to be at least order h p +1 approximations of U and λ∂ U /∂x , respectively, and thus (cid:101) F ∈
Z. Worku, D.W. Zingg P p − ( Ω ) due to the steady version of the primal PDE, (3.1). Similarly, considering (cid:101) ψ, ( (cid:94) λ∂ψ/∂x ) ∈ P p ( Ω ) tobe at least order h p +1 approximations of ψ and ( λ∂ψ/∂x ), respectively, gives (cid:101) G ∈ P p − ( Ω ) due to the adjointPDE, (3.5). For primal and adjoint consistent discretizations, the numerical primal and adjoint solutions are order h ≥ p +1 accurate despite the polynomial approximations; hence, it is sufficient to show that either (4.21) or (4.22)hold for the polynomial integrands instead of the general continuous functions. Note that compatible functionalssatisfy I ( U ) = I ( ψ ), and we enforced the condition I h ( u h ) = I h ( ψ h ) to find the discrete adjoint problem; hence, I ( U ) − I h ( u h ) = I ( ψ ) − I h ( ψ h ). For the rest of the proof, we drop the tilde sign used to distinguish polynomialsfrom the general continuous functions.If U ∈ P p ( Ω ) and ( λ∂ U /∂x ) ∈ P p ( Ω ), then we discretize (3.2) to find I ( U ) = u TL H L g L + g TR H R u R − ψ D w (cid:96)L + ψ N u rR + O (cid:16) h p (cid:17) , (4.23)where w (cid:96)L = [ λ ∂ U ∂x n (cid:96) ] Γ D and u rR = U| Γ N . Subtracting (4.3) from (4.23) and rearranging, we have I ( U ) = I h u h − g TL H L ( u h,L − u L ) + ψ D ( D (cid:96)L u h,L − w (cid:96)L ) − ψ D T ( D ) (cid:96)L ( R (cid:96)L u h,L − u (cid:96)L ) − g TR H R ( u h,R − u R ) − ψ N ( R rR u h,R − u rR ) + O (cid:16) h p (cid:17) . (4.24)Since U ∈ P p ( Ω ), the R γk and D γk matrices are exact when applied to the restriction of U to the grid points, e.g., R rR u R = u rR and D (cid:96)L u L = w (cid:96)L . Applying this property in (4.24) and simplifying we obtain I ( U ) = I h ( u h ) − g T H ( u h − u ) − (cid:20) ψ D T ( D ) (cid:96)L R (cid:96)L − ψ D D (cid:96)L ψ N R rR (cid:21) ( u h − u ) + O (cid:16) h p (cid:17) . (4.25)Adding ψ T H R h,u ( u h ) = 0 to the RHS of (4.25) and rearranging terms, we have I ( U ) = I h ( u h ) − ψ T HD (2) u − ψ T H f + (cid:26) − g T − ψ T HD (2) H − − (cid:20) ψ D T ( D ) (cid:96)L R (cid:96)L − ψ D D (cid:96)L ψ N R rR (cid:21) H − + ψ T AH − + ψ T (cid:34) R T(cid:96)L T ( D ) (cid:96)L R (cid:96)L − D T(cid:96)L R (cid:96)L R TrR D rR (cid:35) H − (cid:27) H ( u h − u ) + O (cid:16) h p (cid:17) . (4.26)Using the identity in (2.6) we can write − HD (2) H − = − (cid:16) D (2) (cid:17) T + (cid:20) D T(cid:96)L R (cid:96)L − R T(cid:96)L D (cid:96)L D TrR R rR − R TrR D rR (cid:21) H − + (cid:20) D TrL R rL − R TrL D rL D T(cid:96)R R (cid:96)R − R T(cid:96)R D (cid:96)R (cid:21) H − + MH − , (4.27)which, after substituting into (4.26) and simplifying, gives I ( U ) = I h ( u h ) − ψ T H (cid:104) D (2) u + f (cid:105) + (cid:26) − g T − ψ T (cid:16) D (2) (cid:17) T + ψ T B T H − + [ F ( ψ )] T H − − ψ T M T H − (cid:27) H ( u h − u ) + O (cid:16) h p (cid:17) . (4.28)Since U ∈ P p ( Ω ), the second term on the RHS vanishes due to the primal PDE. The third term is O ( h ≥ p +1 ) dueto the consistency of the adjoint discretization, the fact that H is O ( h ), and Assumption 1. Therefore, I ( U ) = I h ( u h ) + O (cid:0) h p (cid:1) . Alternatively, if we consider ψ ∈ P p ( Ω ) and ( λ∂ψ/∂x ) ∈ P p ( Ω ), we start by discretizing the second form ofthe functional, (3.4), I ( ψ ) = ψ TL H L f L + ψ TR H R f R − u D z (cid:96)L + u N ψ rR + O (cid:16) h p (cid:17) , (4.29)where z (cid:96)L = (cid:104) λ ∂ψ∂x n (cid:96) (cid:105) Γ D and ψ rR = ψ | Γ N . Subtracting (4.4) from (4.29) and rearranging, we obtain I ( ψ ) = I h ( ψ h ) − f TL H L ( ψ h,L − ψ L ) + u D ( D (cid:96)L ψ h,L − z (cid:96)L ) − u D T ( D ) (cid:96)L ( R (cid:96)L ψ h,L − ψ (cid:96)L ) − f TR H R ( ψ h,R − ψ R ) − u N ( R rR ψ h,R − ψ rR ) + O (cid:16) h p (cid:17) . (4.30)Using the accuracies of R γk and D γk to approximate the boundary terms, adding u T H R h,ψ ( ψ h ) = 0, and simpli-fying leads to I ( ψ ) = I h ( ψ h ) − (cid:104) u T HD (2) ψ + u T H g (cid:105) + (cid:26) − f T − u T HD (2) H − − (cid:20) u D T ( D ) (cid:96)L R (cid:96)L − u D D (cid:96)L u N R rR (cid:21) H − + u T BH − + u T (cid:34) R T(cid:96)L T ( D ) (cid:96)L R (cid:96)L − D T(cid:96)L R (cid:96)L R TrR D rR (cid:35) H − (cid:27) H ( ψ h − ψ ) − u T M ψ h + O (cid:16) h p (cid:17) . (4.31) tability and Functional Superconvergence of Narrow-Stencil SBP Discretizations 9 Using the identity (4.27) and simplifying, we find I ( ψ ) = I h ( ψ h ) − u T H (cid:104) D (2) ψ + g (cid:105) + (cid:26) − f T − u T (cid:16) D (2) (cid:17) T + u T A T H − + [ E ( u )] T H − − u T M T H − (cid:27) H ( ψ h − ψ ) − u T M ψ h + O (cid:16) h p (cid:17) . (4.32)Noting that M + M T = , we have I ( ψ ) = I h ( ψ h ) − u T H (cid:104) D (2) ψ + g (cid:105) + (cid:26) − f T − u T (cid:16) D (2) (cid:17) T + u T A T H − + [ E ( u )] T H − (cid:27) H ( ψ h − ψ ) − u T M ψ + O (cid:16) h p (cid:17) . (4.33)The second and fourth terms on the RHS of (4.33) vanish due to the adjoint PDE and the adjoint consistencyrequirement that M = , respectively. The third term is O ( h ≥ p +1 ) due to the consistency of the primal discretiza-tion, the scaling of the norm matrix, and Assumption 1. Therefore, the estimates in (4.21) and (4.22) hold. (cid:117)(cid:116) Remark 4.1
For implementations with the block-norm wide- or narrow-stencil second-derivative SBP operators ofthe type presented in [27], the estimate in (4.21) is attained even for adjoint inconsistent schemes. Note that forthese types of operator, we have (cid:107) u h,k − u k (cid:107) ∞ = O ( h p ) in (4.28). The block-norm narrow-stencil operators have M k (cid:54) = M Tk ; hence, they lead to adjoint inconsistent schemes even when coupled with adjoint consistent SATs.4.4 Stability AnalysisWe use the energy method to analyze the stability of the SBP-SAT discretization of (3.1). The residual of thediscretization for the homogeneous version of the problem, i.e., F = 0, U D = 0, and U N = 0, summed over allelements can be written as R h ( u h , v ) = − (cid:88) Ω k ∈T h v Tk M k u h,k − (cid:88) γ ⊂ Γ D (cid:20) R γk v k D γk v k (cid:21) T (cid:34) T ( D ) γk − − (cid:35) (cid:20) R γk u h,k D γk u h,k (cid:21) − (cid:88) γ ⊂ Γ I R γk v k R γv v v D γk v k D γv v v T T (1) γk − T (1) γk T (3) γk − T (3) γk − T (1) γv T (1) γv T (3) γv T (3) γv − T (2) γk − T (2) γk T (4) γk T (4) γk − T (2) γv T (2) γv T (4) γv T (4) γv R γk u h,k R γv u h,v D γk u h,k D γv u h,v (4.34)for v ∈ R n e n p . In [37,36], a factorization of M k for wide-stencil operators allowed the use of the borrowing trickand enabled R h ( u h , v ) to be written in terms of interface contributions only. However, the same factorization cannot be applied for narrow-stencil operators because D γk is constructed using a modified derivative operator atthe element boundaries, D b,k , instead of D k . To circumvent this, we make the following assumption. Assumption 2
The D b,k matrix is invertible or can be modified such that it is invertible. The invertibility requirement on D b,k is not too restrictive. In fact, all the narrow-stencil second-derivativeoperators in [29,28,26,27] either have invertible D b,k matrix or their D b,k matrix can be modified such that it isinvertible. For operators that include nodes at element boundaries, the only requirement for D b,k to be invertible isthat its interior diagonal entries are nonzero, e.g., D b,k can be constructed from the identity matrix by modifyingthe first and last rows such that these rows approximate the first derivative to degree ≥ p + 1. The invertibility of D b,k matrix constructed in this manner can easily be verified using Gershgorin’s theorem. For generalized narrow-stencil second-derivative operators with nodes at element boundaries, e.g., the hybrid Gauss-trapezoidal-Lobatto(HGTL) operators in [12], a similar modification can be applied to obtain an invertible D b,k matrix. In contrast,all except the degree two hybrid Gauss-trapezoidal (HGT) operators in [12], which do not include boundary nodes,do not yield an invertible D b,k matrix even after applying the modification discussed. However, it is likely possibleto construct HGT operators such that D b,k is invertible by enforcing a condition on the free variables during the construction of the operators. The structures of the invertible D b,k matrices of the degree two CSBP and HGT operators are, × × × × × × × × , × × × × ×× × × × ×× × × × ×× × × × × × × × × ×× × × × ×× × × × ×× × × × × , respectively, where each row containing × in its entries approximate the first derivative.Another important assumption that is required in the subsequent energy stability analysis for adjoint consistentdiscretizations with narrow-stencil second-derivative operators is presented below. Assumption 3
The first and second-derivative SBP operators, D k and D (2) k , are nullspace consistent, i.e., thenonzero vectors in the nullspace of D k and D (2) k are N ( D k ) = span { } =: v c and N ( D (2) k ) = span { , x k } , respec-tively. It should be noted that consistency of an SBP operator does not necessarily imply nullspace consistency and viceversa. The operators defined in Definitions 2.1 and 2.2 are consistent because they satisfy the accuracy conditions,i.e., they differentiate polynomials up to a required degree exactly [35]. In contrast, nullspace consistency requiresthat the nullspaces of D k and D (2) k exclusively contain vectors in span { } and span { , x k } , respectively. SBPderivative operators are consistent by construction, and most of them are nullspace consistent as well [35].Using Assumption 2 and enforcing the conditions necessary for conservation, (4.14), we can now write the sumof the residual and its transpose as2 R h ( u h , v h ) = R h ( u h , v h ) + R Th ( u h , v h ) = − (cid:88) γ ⊂ Γ I R γk u h,k R γv u h,v D b,k u h,k D b,v u h,v T T (1) γk − T (1) γk σ k C γk − σ v C γv − T (1) γk T (1) γv − σ k C γk σ v C γv σ k C Tγk − σ k C Tγk α γk V k − σ v C Tγv σ v C Tγv α γv V v R γk u h,k R γv u h,v D b,k u h,k D b,v u h,v − (cid:88) γ ⊂ Γ I (cid:20) D γk u h,k D γv u h,v (cid:21) T (cid:34) T (4) γk T (4) γk T (4) γk T (4) γ (cid:35) (cid:20) D γk u h,k D γv u h,v (cid:21) − (cid:88) γ ⊂ Γ D (cid:20) R γk u h,k D b,k u h,k (cid:21) T (cid:34) T ( D ) γk − C γk − C Tγk α γk V k (cid:35) (cid:20) R γk u h,k D b,k u h,k (cid:21) , (4.35)where C γk = n γk R γk Λ k , C γv = n γv R γv Λ v , σ k = T (2) γk + T (3) γk − σ v = T (2) γv + T (3) γv − α γk is a positive interfaceweight factor satisfying the relation (cid:80) γ ∈ Γ k α γk = 1, and V k = D − Tb,k ( M k + M Tk ) D − b,k , V v = D − Tb,v ( M v + M Tv ) D − b,v . (4.36)We note that V k is positive semidefinite since v T ( M k + M Tk ) v ≥ v ∈ R n p implies( D − b,k v ) T ( M k + M Tk )( D − b,k v ) ≥ . (4.37)Moreover, we have D b,k v c = v , or D − b,k v = v c , (4.38)where v represents vectors containing zero at the entries corresponding to the rows for which D b,k containsconsistent approximations of the first derivative and the values of v c at all other entries.For diagonal-norm narrow-stencil SBP operators that are constructed as in [29,28,26,12,27], we can determinethe vectors in the nullspace of V k using Assumptions 2 and 3. Lemma 4.1
Consider a consistent diagonal-norm narrow-stencil second-derivative SBP operator of the form (2.4) for which M k = M Tk , the E k matrix is constructed such that it has nonzero rows only at row indices where the D b,k matrix contains consistent approximations of the first derivative, the D (2) k matrix has larger dense blocks at the topleft and bottom right corners than the E k matrix, and Assumptions 2 and 3 hold. Then, we have N ( M k ) = v c and N ( V k ) = N ( V v ) = v . tability and Functional Superconvergence of Narrow-Stencil SBP Discretizations 11 Proof
The nullspace consistency of the D (2) k in Assumption 3 implies the following: H k D (2) k v c = ( − M k + E k Λ k D b,k ) v c = − M k v c + E k Λ k D b,k v c = , (4.39) H k D (2) k x k = ( − M k + E k Λ k D b,k ) x k = − M k x k + E k Λ k D b,k x k = . (4.40)The second term in the last equality in (4.39) is zero due to the structure of the E k matrix, i.e., E k Λ k D b,k v c = ;thus, M k v c = . Furthermore, M k x k (cid:54) = in (4.40) since otherwise we would obtain H k D (2) k x k = E k Λ k D b,k x k = R Trk D rk x k + R T(cid:96)k D (cid:96)k x k = R Trk λ r − R T(cid:96)k λ (cid:96) = , (4.41)which is not possible as R rk and R (cid:96)k do not have nonzero values at the same entries, and λ >
0. We have used theaccuracy of D rk and D (cid:96)k in the penultimate equality in (4.41), i.e., D rk x k = λ r and D (cid:96)k x k = − λ (cid:96) , where λ (cid:96) and λ r are at least order h p +1 approximations of λ at the left and right boundaries of Ω k , respectively. Hence, thereis no vector spanned by { , x k } other than v c that is in the nullspace of M k . If there exists a nontrivial vector v such that v / ∈ span { , x k } and M k v = , then H k D (2) k v = E k Λ k D b,k v (cid:54) = , (4.42)because D (2) k is nullspace consistent and H k is SPD. The vector E k Λ k D b,k v has zero entries at rows correspondingto the zero rows of the E k matrix. By construction, D (2) k has larger dense blocks at the top left and bottomright corners (consisting of more rows and columns) than the E k matrix; therefore, it follows from the nullspaceconsistency of the D (2) k matrix that [ H k D (2) k v ] i (cid:54) = 0 and [ E k Λ k D b,k v ] i = 0, at least for one entry, the i -th entry,near the boundaries. This implies that the equality in (4.42) cannot hold for any vector v / ∈ span { , x k } ; hence,we have N ( M k ) = v c . Since M k = M Tk , it follows that N ( M Tk ) = v c . Using the result in (4.38) with the fact that N ( M k ) = N ( M Tk ) = v c , we obtain V k v = D − Tb,k ( M k + M Tk ) D − b,k v = D − b,k ( M k + M Tk ) v c = . (4.43)Thus, v is the only nontrivial vector in the nullspace of V k . Analogous results hold for V v . (cid:117)(cid:116) In [14], the stability conditions that the SATs must satisfy were derived for diagonal-norm narrow-stencilSBP operators assuming that M k is SPD; however, most operators in the literature, e.g., [29,28,26,12,27], do notsatisfy this requirement. For dense-norm narrow-stencil second-derivative SBP operators, we make the followingassumption regarding the nullspaces of M k and M Tk : Assumption 4
For dense-norm narrow-stencil SBP operators, v c is the only nontrivial vector in the nullspacesof M k and M Tk , i.e., N ( M k ) = N ( M Tk ) = v c . Under Assumption 4, (4.43) gives N ( V k ) = N ( V v ) = v for dense-norm narrow-stencil SBP operators. Beforeproceeding with the energy analysis of the SBP-SAT discretization, we state an essential theorem, which is provedin [1,15]. Theorem 4.2
A symmetric matrix of the form Y = (cid:2) Y Y Y T Y (cid:3) is positive semidefinite if and only if Y (cid:23) , ( I − Y Y +22 ) Y T = , and Y − Y Y +22 Y T (cid:23) , (4.44) where Y + denotes the Moore-Penrose pseudoinverse of Y and Y (cid:23) indicates that Y is positive semidefinite. An SBP-SAT discretization is energy stable ifdd t (cid:16) (cid:107) u h (cid:107) H (cid:17) = u Th H d u h d t + d u Th d t H u h = 2 R h ( u h , u h ) ≤ . (4.45)The sum of the residual and its transpose for conservative schemes, 2 R h ( u h , u h ), given in (4.35) satisfies the energystability condition if A := (cid:20) A A A A (cid:21) = T (1) γk − T (1) γk σ k C γk − σ v C γv − T (1) γk T (1) γv − σ k C γk σ v C γv σ k C Tγk − σ k C Tγk α γk V k − σ v C Tγv σ v C Tγv α γv V v , (cid:34) T (4) γk T (4) γk T (4) γk T (4) γ (cid:35) , (cid:34) T ( D ) γk − C γk − C Tγk α γk V k (cid:35) , (4.46)are positive semidefinite. We partition the matrix A ∈ R (2+2 n p ) × (2+2 n p ) using four blocks, namely A ∈ R × , A ∈ R × n p , A ∈ R n p × , and A ∈ R n p × n p .For adjoint consistent SATs, we enforce all the conditions in (4.11). Furthermore, to obtain a symmetric A matrix we require that T (3) γk − T (2) γk = 1, as in [37]. This condition is satisfied by the SATs corresponding to some of the popular discontinuous Galerkin fluxes for elliptic PDEs [36], e.g., the modified method of Bassi and Rebay (BR2) [3], local discontinuous Galerkin (LDG) [33], and compact discontinuous Galerkin (CDG) [31] methods.With the adjoint consistency and symmetry conditions in place, the components of the A matrix become A = (cid:34) T (1) γk − T (1) γk − T (1) γk T (1) γv (cid:35) , A = (cid:34) T (2) γk C γk − T (2) γv C γv − T (2) γk C γk T (2) γv C γv (cid:35) , A = (cid:34) T (2) γk C Tγk − T (2) γk C Tγk − T (2) γv C Tγv T (2) γv C Tγv (cid:35) , A = (cid:20) α γk V k α γv V v (cid:21) . (4.47) Theorem 4.3
Let Assumptions 2 to 4 hold, then an adjoint consistent SBP-SAT discretization of the diffusionproblem, (3.1) , that produces a symmetric A matrix, i.e., T (3) γk − T (2) γk = 1 , and uses a narrow-stencil second-derivativeoperator of the form (2.4) for the spatial discretization is energy stable if T (1) γk ≥ α γk T (2) γk R γk Λ k V + k Λ k R Tγk T (2) γk + 2 α γv T (2) γv R γv Λ v V + v Λ v R Tγv T (2) γv , (4.48) T ( D ) γk ≥ α γk R γk Λ k V + k Λ k R Tγk , (4.49) T (4) γk ≥ . (4.50) Proof
We wish to show that the matrices in (4.46) are positive semidefinite. The matrix A , whose componentsare given in (4.47), is symmetric; thus, we can use Theorem 4.2 to determine the conditions required for it tobe positive semidefinite. We have A (cid:23) V k and V v are positive semidefinite, α γk >
0, and α γv > (cid:16) I − A A +22 (cid:17) A = 2 (cid:20) I − V k V + k I − V v V + v (cid:21) (cid:34) C Tγk T (2) γk − C Tγk T (2) γk − C Tγv T (2) γv C Tγv T (2) γv (cid:35) = 0 . (4.51)To show that (4.51) holds, we consider the singular value decomposition of V k , V k = X Σ Y T , (4.52)where the columns of X and Y contain orthonormal basis vectors of the column and row spaces, respectively, and Σ is a diagonal matrix containing the singular values of V k along its diagonal. Lemma 4.1 and Assumption 4 ensurethat the matrix V k ∈ R n p × n p has only one nontrivial vector in its nullspace; hence, the first n p − X contain orthonormal basis vectors that span the column space of V k and the last column contains the vector inthe nullspace of V k , which is v . We also note that V k V + k = X Σ Y T Y Σ + X T = X ΣΣ + X T = XI m X T , (4.53)where we have used the orthonormality of Y in the second equality, i.e., Y T Y = I , and I m denotes the identitymatrix with the last diagonal entry set to zero (since Σ ii = 0 for i > m , where m = n p − V k ).Therefore, for operators that include the boundary nodes, we have I − V k V + k = I − XI m X T , which gives I − XI m X T = I − × . . . × × . . . × × .. . ... .. . × . . . × ×× . . . × . . . 1 0 × × . . . × × ... ... ... ... × × . . . × × × . . . × = . . . × . . . × . ... .. . ...0 × . . . ×
00 0 . . . , (4.54)where × denotes an entry that we do not need to specify for this analysis. Similarly, for operators that do notinclude boundary nodes, we obtain I − V k V + k = I − XI m X T = × . . . × . .. . .. × . . . × , (4.55) tability and Functional Superconvergence of Narrow-Stencil SBP Discretizations 13 i.e., the first and last s rows are zero, where s is half of the number of zero entries in v . For operators that havenodes at the boundaries, the LHS of (4.51) can be evaluated using (4.54) as (cid:16) I − A A +22 (cid:17) A = . . . × . . . × × . . . ×
00 0 . . . . . . × . . . × × . . . ×
00 0 . . . T (2) γk λ γk − T (2) γk λ γk − T (2) γv λ γv T (2) γv λ γv = 0 . (4.56)Similarly, substituting (4.55) into (4.51), it is straightforward to show that the condition (cid:0) I − A A +22 (cid:1) A = 0also holds for operators that do not include boundary nodes. For A to be positive semidefinite, it remains to findsufficient conditions to satisfy the last requirement in Theorem 4.3, i.e., A − A A +22 A (cid:23)
0, which, after somealgebra, gives the condition (cid:20) − − (cid:21) ⊗ (cid:20) T − (cid:18) α γk T (2) γk C γk V + k C Tγk T (2) γk + 2 α γv T (2) γv C γv V + v C Tγv T (2) γv (cid:19)(cid:21) (cid:23) , (4.57)where ⊗ denotes the Kronecker product. Since (cid:2) − − (cid:3) (cid:23)
0, the inequality in (4.57) is satisfied if T (1) γ ≥ α γk T (2) γk C γk V + k C Tγk T (2) γk + 2 α γv T (2) γv C γv V + v C Tγv T (2) γv , (4.58)which is the same as the condition given by (4.48). Note that C γk V + k C Tγk = R γk Λ k V + k Λ k R Tγk since n γk = 1.The second matrix in (4.46) can be written as (cid:34) T (4) γk T (4) γk T (4) γk T (4) γk (cid:35) = 2 T (4) γk (cid:20) (cid:21) , (4.59)which is positive semidefinite provided T (4) γk = T (4) γv ≥
0, since (cid:2) (cid:3) (cid:23) α γk V k (cid:23)
0. Furthermore, using the same approach used to obtain (4.56) it can be shownthat ( I − V k V + k )( − C Tγk ) = 0 (4.60)irrespective of whether or not the operator includes boundary nodes. The last condition required for positivesemidefiniteness of the last matrix in (4.46) is T ( D ) γk − C γk ( α γk V k ) + C Tγk ≥ , (4.61)which is satisfied if (4.49) holds. Therefore, the conditions in Theorem 4.3 are indeed sufficient for energy stability. (cid:117)(cid:116) Remark 4.2
In practice, the pseudoinverse V + k for a given SBP operator is computed once and for all on thereference element. It is then scaled by the inverse of the metric Jacobian when the operator is mapped to thephysical elements. In Table 1 of [13], equivalent values of 2 R γk V + k R Tγk , denoted by q and scaled by the meshspacing, are tabulated for the constant-coefficient diagonal-norm narrow-stencil SBP operators presented in [29].The scaling used can be written as h = 1 / [( n e n p − − ( n e − V + k = ( D − Tb,k M k D − b,k ) + for operators with M k = M Tk is computed as D b,k ( (cid:101) M k ) − D Tb,k , where (cid:101) M k is obtained by perturbing M k such that the corner values of D b,k ( (cid:101) M k ) − D Tb,k are independent of the perturbation. The difference in the values of qh and 2 h R γk V + k R Tγk lies inthe approaches pursued to evaluate 2 V + k . Remark 4.3
For stability of discretizations with wide-stencil second-derivative operators, the terms V + k and V + v in Theorem 4.3 are replaced by ( H k Λ k + Λ k H k ) − and ( H v Λ v + Λ v H v ) − , respectively. A stabilized version of the BR2 method [3] for implementation with the narrow-stencil second-derivative SBPoperators can be obtained by choosing T (1) γk = T (1) γv = 12 α γk R γk Λ k V + k Λ k R Tγk + 12 α γv R γv Λ v V + v Λ v R Tγv , − T (2) γk = − T (2) γv = T (3) γk = T (3) γk = 12 , T ( D ) γk = 2 α γk R γk Λ k V + k Λ k R Tγk , T (4) γk = T (4) γv = 0 . (4.62)The general form of the BR2 SAT was first proposed in [37] by discretizing the primal formulation of the DGmethod using SBP operators. It is straightforward to show that the BR2 SAT coefficients satisfy the adjointconsistency conditions in (4.11) and the stability requirements in Theorem 4.3. We determine the SAT coefficients corresponding to the LDG method [33] by discretizing the primal LDG for-mulation of the diffusion problem (see, e.g., [2,31] for the primal LDG formulation). Similar analysis with thewide-stencil second-derivative SBP operators can be found in [8,16,5,36]. Stable LDG SAT coefficients (with nomesh dependent parameter for stabilization) for discretizations with the narrow-stencil second-derivative SBPoperators are given by T (1) γk = T (1) γv = T ( D ) γk = 2 α γk R γk Λ k V + k Λ k R Tγk , − T (2) γk = T (3) γv = 1 , T (3) γk = T (2) γv = T (4) γk = T (4) γv = 0 . (4.63)Clearly, the LDG SAT coefficients in (4.63) satisfy both the adjoint consistency conditions in (4.11) and thestability demands in Theorem 4.3. Remark 4.4
The LDG SAT coefficients presented in (4.63) are obtained by using a switch function value of 1 / x -axis. We refer the reader to [2,31,36] for details regarding the switchfunction and to [32] for a discussion on the need to use a global vector. Remark 4.5
In one space dimension the LDG and CDG fluxes are identical [31]; therefore, the SAT coefficientspresented in (4.63) define the CDG SAT as well.
The SAT coefficients corresponding to the BO method [4] do not satisfy the adjoint consistency conditions in (4.11);hence, the energy stability requirements in Theorem 4.3 do not apply. The BO SAT coefficients for implementationwith narrow-stencil second-derivative operators are given by T (2) γk = T (2) γv = T (3) γk = T (3) γv = 12 , T ( D ) γk = 2 α γk R γk Λ k V + k Λ k R Tγk , T (1) γk = T (1) γv = T (4) γk = T (4) γv = 0 . (4.64)Except for the conditions on T ( D ) γk , the stability conditions for the narrow- and wide-stencil SBP operators are thesame when the BO SAT is used. Stability analysis for discretizations with wide-stencil SBP operators and the BOSAT can be found in [8,16,36]. Note that for the BO SAT coefficients in (4.64), we have T (1) γk = T (1) γv = σ k = σ v = 0in the matrix A given in (4.46), and thus A is positive semidefinite. The stability analyses for the second and thirdmatrices in (4.46) remain the same as those presented in the proof of Theorem 4.3. The BO SAT satisfies the conditions for conservation (4.14); hence, it leads to a conservative and stable but not adjoint consistent scheme. tability and Functional Superconvergence of Narrow-Stencil SBP Discretizations 15 A version of the CNG SAT [7] that leads to a stable discretization when implemented with narrow-stencil second-derivative SBP operators has the SAT coefficients T (1) γk = T (1) γv = 18 α γk R γk Λ k V + k Λ k R Tγk + 18 α γv R γv Λ v V + v Λ v R Tγv , T (2) γk = T (2) γv = T (4) γk = T (4) γv = 0 , T (3) γk = T (3) γv = 12 , T ( D ) γk = 2 α γk R γk Λ k V + k Λ k R Tγk . (4.65)Clearly, the coefficients in (4.65) satisfy all the conditions in (4.11) except the second one. Therefore, the CNGSAT leads to conservative but adjoint inconsistent schemes. The stability analyses of the second and third matricesin (4.46) are the same as those presented in the proof of Theorem 4.3. The positive semidefiniteness of the matrix A in (4.46) requires all the conditions in (4.44) to be satisfied. Substituting the CNG SAT coefficients in A , we seethat A (cid:23)
0, and it can be shown that (cid:0) I − A A +22 (cid:1) A = 0. Hence, it only remains to find conditions such that A − A A +22 A (cid:23)
0, which, after simplification, yields (cid:20) − − (cid:21) ⊗ (cid:20) T (1) γk − (cid:18) α γk C γk V + k C Tγk + 14 α γk C γv V + v C Tγv (cid:19)(cid:21) (cid:23) . (4.66)The inequality in (4.66) is satisfied by the T (1) γk = T (1) γv coefficient given in (4.65). We consider the one-dimensional Poisson problem with Dirichlet and Neumann boundary conditions, − ∂ U ∂x = F in Ω = [0 , , U (cid:12)(cid:12) x =0 = U D , ∂ U ∂x (cid:12)(cid:12)(cid:12)(cid:12) x =1 = U N . (5.1)We use the method of manufactured solution and let U = cos(30 x ) as in [13]; thus, F = 30 cos(30 x ). We alsoconsider a compatible linear functional given by I ( U ) = (cid:90) cos (30 x ) d Ω + 130 (1 −
30 sin(30) − cos(30)) cos(30) . (5.2)We are interested in the convergence of the solution and functional errors under mesh refinement. Figure 5.1presents the solution convergence for discretizations with the diagonal-norm narrow-stencil CSBP operators in [29]and the generalized SBP operators in [12] that have an invertible D b,k matrix and satisfy the accuracy conditions,i.e., the p = { , } HGTL and p = 2 HGT operators. Note that the degree four HGTL operator in [12] meetsthe accuracy requirements given in Definitions 2.1 and 2.2 to order h only, while the degree three and four HGToperators have D b,k matrices that cannot be modified as described in Section 4 to ensure their invertibility. Theinterface weight parameters in the SAT coefficients are set as α γk = α γv = 1 / The solution error is computed as (cid:115) (cid:88) Ω k ∈T h ( u h,k − u k ) T H k ( u h,k − u k ) . The convergence rates in Fig. 5.1 through Fig. 5.8 are calculated by fitting a line through the error values on themesh resolutions indicated by the short, thin lines, and “dof” stands for the number of degrees of freedom in thespatial discretization. Figure 5.1 shows that a solution convergence rate of p + 2 is attained when order-matchednarrow-stencil operators are coupled with the adjoint consistent SATs, except with the degree one SBP operators.The adjoint inconsistent SATs, BO and CNG, exhibit solution convergence rates of p+2 with all the order-matchednarrow-stencil SBP operators, except the degree one and three CSBP operators which yield convergence rates of p+ 1. This is consistent with the results presented in [12] but somewhat surprising since the well-known even-oddconvergence phenomenon (see e.g., [8,33,25]) that the BO method displays is not observed with the narrow-stencilSBP operators. Numerical experiments in [36] show that the BO and CNG SATs converge at rates of p + 1 and p when implemented with odd and even degree multidimensional SBP operators, respectively. Indeed, this even-odd Although not presented here, we observe a solution convergence rate of p + 2 with all the diagonal-norm narrow-stencil CSBPoperators presented in [12].6 Z. Worku, D.W. Zingg(a) p = 1 (b) p = 2(c) p = 3 (d) p = 4 Fig. 5.1
Solution convergence under mesh refinement. The values in parentheses are the convergence rates, and “N” stands fornarrow-stencil SBP operator. convergence phenomenon is also observed when the BO and CNG SATs are implemented with the Legendre-Gauss-Lobatto (LGL) and Legendre-Gauss (LG) wide-stencil SBP operators, as shown in Fig. 5.2. However, thistrend does not hold consistently with the wide-stencil CSBP and HGT operators, as convergence rates of p + 1 areachieved with the p = 4 operators, as depicted in Fig. 5.2d. Therefore, it appears that the even-odd convergenceproperty of the BO and CNG SATs is dependent on the type of SBP operator used, and it is not observed withthe diagonal-norm narrow-stencil SBP operators consistently. Implementations of the diagonal-norm wide-stencilSBP operators with the BR2 and LDG SATs lead to solution convergence rates of p + 1, as depicted in Fig. 5.3.Finally, Fig. 5.4 shows that the block-norm wide- and narrow-stencil SBP operators presented in [27], denoted byCSBP2, achieve a solution convergence rate of 2 p regardless of the type of SAT used.The functional error is calculated as | I h ( u h ) −I ( U ) | . Figure 5.5 shows the functional convergence rates resultingfrom discretizations with the order-matched narrow-stencil SBP operators. As established in Theorem 4.1, thefigure shows that the functional superconverges at a rate of 2 p when adjoint consistent SATs are used. The adjoint inconsistent SATs yield larger functional error values and lower functional convergence rates when coupledwith the degree three and four diagonal-norm narrow-stencil SBP operators, but they attain lower error valuesand a convergence rate of 2 p for degree one and two operators. As can be seen from Fig. 5.6, when the adjointinconsistent SATs are used with the diagonal-norm wide-stencil SBP operators, convergence rates of 2 p are notattained, except for the p = 1 case. It is also evident from Figs. 5.5 and 5.6 that in most cases the diagonal-normnarrow-stencil operators result in lower functional error and larger functional convergence rates than the diagonal-norm wide-stencil operators when used with the adjoint inconsistent SATs. Figure 5.7 shows that the adjointconsistent SATs lead to functional convergence rates of 2 p when used with the diagonal-norm wide-stencil SBPoperators. Comparing the results depicted in Figs. 5.5 and 5.7, we can conclude that diagonal-norm narrow-stencilSBP operators do not offer better functional convergence rates than diagonal-norm wide-stencil SBP operatorswhen coupled with the adjoint consistent SATs, which agrees with the theory. Similarly, the functional convergencerates attained with the block-norm wide- and narrow-stencil SBP operators are comparable, as depicted in Fig. 5.8.Furthermore, the functional convergence rate with the block-norm SBP operators is 2 p for adjoint consistent as wellas adjoint inconsistent schemes, except for the degree three block-norm wide-stencil SBP operator, which exhibitsa 2 p − narrow-stencil SBP operators. tability and Functional Superconvergence of Narrow-Stencil SBP Discretizations 17(a) p = 1 (b) p = 2(c) p = 3 (d) p = 4 Fig. 5.2
Solution convergence under mesh refinement with adjoint inconsistent SATs. The values in parentheses are the convergencerates, and “W” stands for wide-stencil SBP operator.(a) p = 1 (b) p = 2(c) p = 3 (d) p = 4 Fig. 5.3
Solution convergence under mesh refinement with adjoint consistent SATs. The values in parentheses are the convergencerates, and “W” stands for wide-stencil SBP operator.8 Z. Worku, D.W. Zingg(a) p = 1 (b) p = 2(c) p = 3 (d) p = 4 Fig. 5.4
Solution convergence under mesh refinement with block-norm wide-stencil (“W”) and narrow-stencil (“N”) CSBPoperators. The values in parentheses are the convergence rates.(a) p = 1 (b) p = 2(c) p = 3 (d) p = 4 Fig. 5.5
Functional convergence under mesh refinement. The values in parentheses are the convergence rates, and “N” stands fornarrow-stencil SBP operator.tability and Functional Superconvergence of Narrow-Stencil SBP Discretizations 19(a) p = 1 (b) p = 2(c) p = 3 (d) p = 4 Fig. 5.6
Functional convergence under mesh refinement with adjoint inconsistent SATs. The values in parentheses are theconvergence rates, and “W” stands for wide-stencil SBP operator. The convergence rates for the p = 3 HGTL and p = 4 HGToperators are calculated by fitting lines through the error values on the mesh resolution indicated by the short, thin lines with astar marker. In this paper, we have shown that primal and adjoint consistent SBP-SAT discretizations of diffusion problemswith diagonal-norm second-derivative generalized SBP operators lead to functional superconvergence if the primaland adjoint solutions are sufficiently smooth. For block-norm second-derivative operators, however, the analysisand the numerical experiments show that adjoint inconsistency does not degrade the functional convergence rate.We have also derived the conditions required for the stability of adjoint consistent SBP-SAT discretizations withthe narrow-stencil second-derivative generalized SBP operators under the assumptions that the operators areconsistent and nullspace consistent. The stability analysis also requires that the derivative operator at the elementboundaries, D b,k , be invertible. For most operators, D b,k is invertible or can easily be modified to be invertible. Forsome operators, however, this is not the case, and it might be necessary to enforce the invertibility of this matrix during the construction of the SBP operators to ensure that SBP-SAT discretizations with these operators arestable and adjoint consistent in addition to the other attractive numerical properties that narrow-stencil operatorsoffer.Four different types of stable SATs for narrow-stencil SBP operators, among which two are adjoint consistent,are proposed and implemented in the numerical experiments. As predicted by the theory, the numerical experimentsshow that functionals superconverge at a rate of 2 p when a diagonal-norm degree p + 1 narrow-stencil or degree p wide-stencil generalized SBP operator is used along with adjoint consistent SATs. It is also observed that theadjoint consistent BR2 and LDG SATs yield solution convergence rates of p + 1 and p + 2 when implemented withdiagonal-norm wide- and narrow-stencil SBP operators, respectively. Implementations with the block-norm wide-and narrow-stencil SBP operators show a solution and functional convergence rates of 2 p regardless of the type ofSAT used. The even-odd convergence properties of the adjoint inconsistent BO and CNG SATs are not observedwhen these SATs are implemented with the diagonal- and block-norm narrow-stencil SBP operators.While the SATs presented in this work ensure the consistency, conservation, adjoint consistency, stability, andfunctional superconvergence of SBP-SAT discretizations with narrow-stencil generalized SBP operators, optimiza-tion of the SAT coefficients to achieve improved numerical properties, e.g., better conditioning and spectral radius, may be pursued in future work. p = 1 (b) p = 2(c) p = 3 (d) p = 4 Fig. 5.7
Functional convergence under mesh refinement with adjoint consistent SATs. The values in parentheses are the conver-gence rates, and “W” stands for wide-stencil SBP operator.(a) p = 1 (b) p = 2(c) p = 3 (d) p = 4 Fig. 5.8
Functional convergence under mesh refinement with block-norm wide-stencil (“W”) and narrow-stencil (“N”) CSBPoperators. The values in parentheses are the convergence rates.tability and Functional Superconvergence of Narrow-Stencil SBP Discretizations 21
References
1. Albert, A.: Conditions for positive and nonnegative definiteness in terms of pseudoinverses. SIAM Journal on Applied Mathe-matics (2), 434–440 (1969)2. Arnold, D.N., Brezzi, F., Cockburn, B., Marini, L.D.: Unified analysis of discontinuous Galerkin methods for elliptic problems.SIAM Journal on Numerical Analysis (5), 1749–1779 (2002)3. Bassi, F., Rebay, S., Mariotti, G., Pedinotti, S., Savini, M.: A high-order accurate discontinuous finite element method forinviscid and viscous turbomachinery flows. In: R. Decuypere, G. Dibelius (eds.) Proceedings of the 2nd European Conferenceon Turbomachinery, Fluid Dynamics and Thermodynamics, pp. 99–109. Technologisch Instituut, Antwerpen, Belgium (1997)4. Baumann, C.E., Oden, J.T.: A discontinuous hp finite element method for convection–diffusion problems. Computer Methodsin Applied Mechanics and Engineering (3-4), 311–341 (1999)5. Berg, J., Nordstr¨om, J.: Superconvergent functional output for time-dependent problems using finite differences on Summation-By-Parts form. Journal of Computational Physics (20), 6846–6860 (2012)6. Carpenter, M.H., Gottlieb, D., Abarbanel, S.: Time-stable boundary conditions for finite-difference schemes solving hyperbolicsystems: methodology and application to high-order compact schemes. Journal of Computational Physics (2), 220–236(1994)7. Carpenter, M.H., Nordstr¨om, J., Gottlieb, D.: A stable and conservative interface treatment of arbitrary spatial accuracy.Journal of Computational Physics (2), 341–365 (1999)8. Carpenter, M.H., Nordstr¨om, J., Gottlieb, D.: Revisiting and extending interface penalties for multi-domain summation-by-parts operators. Journal of Scientific Computing (1-3), 118–150 (2010)9. Craig Penner, D.A., Zingg, D.W.: Superconvergent functional estimates from tensor-product generalized summation-by-partsdiscretizations in curvilinear coordinates. Journal of Scientific Computing (2), 41 (2020)10. Del Rey Fern´andez, D.C., Boom, P.D., Zingg, D.W.: A generalized framework for nodal first derivative summation-by-partsoperators. Journal of Computational Physics , 214–239 (2014)11. Del Rey Fern´andez, D.C., Hicken, J.E., Zingg, D.W.: Simultaneous approximation terms for multi-dimensional summation-by-parts operators. Journal of Scientific Computing (1), 83–110 (2018)12. Del Rey Fern´andez, D.C., Zingg, D.W.: Generalized summation-by-parts operators for the second derivative. SIAM Journalon Scientific Computing (6), A2840–A2864 (2015)13. Eriksson, S.: A dual consistent finite difference method with narrow stencil second derivative operators. Journal of ScientificComputing (2), 906–940 (2018)14. Eriksson, S., Nordstr¨om, J.: Finite difference schemes with transferable interfaces for parabolic problems. Journal of Compu-tational Physics , 935–949 (2018)15. Gallier, J.H.: Notes on the Schur complement. Penn Engineering pp. 1–12 (2010)16. Gong, J., Nordstr¨om, J.: Interface procedures for finite difference approximations of the advection–diffusion equation. Journalof Computational and Applied Mathematics (5), 602–620 (2011)17. Gustafsson, B.: The convergence rate for difference approximations to general mixed initial-boundary value problems. SIAMJournal on Numerical Analysis (2), 179–190 (1981)18. Hartmann, R.: Adjoint consistency analysis of discontinuous Galerkin discretizations. SIAM Journal on Numerical Analysis (6), 2671–2696 (2007)19. Hartmann, R., Leicht, T.: Higher order and adaptive DG methods for compressible flows. In: Deconinck H, ed. VKI LS 2014-03: 37th Advanced VKI CFD Lecture Series: Recent Developments in Higher Order Methods and Industrial Application inAeronautics, Dec. 9-12, 2013. Rhode-Saint-Genese, Belgium: Von Karman Institute for Fluid Dynamics (2014)20. Hicken, J.E.: Output error estimation for summation-by-parts finite-difference schemes. Journal of Computational Physics (9), 3828–3848 (2012)21. Hicken, J.E., Del Rey Fern´andez, D.C., Zingg, D.W.: Multidimensional summation-by-parts operators: General theory andapplication to simplex elements. SIAM Journal on Scientific Computing (4), A1935–A1958 (2016)22. Hicken, J.E., Zingg, D.W.: Superconvergent functional estimates from summation-by-parts finite-difference discretizations.SIAM Journal on Scientific Computing (2), 893–922 (2011)23. Hicken, J.E., Zingg, D.W.: Summation-by-parts operators and high-order quadrature. Journal of Computational and AppliedMathematics (1), 111–125 (2013)24. Hicken, J.E., Zingg, D.W.: Dual consistency and functional accuracy: a finite-difference perspective. Journal of ComputationalPhysics , 161–182 (2014)25. Kirby, R.M., Karniadakis, G.E.: Selecting the numerical flux in discontinuous Galerkin methods for diffusion problems. Journalof Scientific Computing (1-3), 385–411 (2005)26. Mattsson, K.: Summation by parts operators for finite difference approximations of second-derivatives with variable coefficients.Journal of Scientific Computing (3), 650–682 (2012)27. Mattsson, K., Almquist, M.: A solution to the stability issues with block norm summation by parts operators. Journal ofComputational Physics , 418–442 (2013)28. Mattsson, K., Ham, F., Iaccarino, G.: Stable and accurate wave-propagation in discontinuous media. Journal of ComputationalPhysics (19), 8753–8767 (2008)29. Mattsson, K., Nordstr¨om, J.: Summation by parts operators for finite difference approximations of second derivatives. Journalof Computational Physics (2), 503–540 (2004)30. Mattsson, K., Sv¨ard, M., Shoeybi, M.: Stable and accurate schemes for the compressible Navier–Stokes equations. Journal ofComputational Physics (4), 2293–2316 (2008)31. Peraire, J., Persson, P.O.: The compact discontinuous Galerkin (CDG) method for elliptic problems. SIAM Journal on ScientificComputing (4), 1806–1824 (2008)32. Sherwin, S., Kirby, R., Peir´o, J., Taylor, R., Zienkiewicz, O.: On 2D elliptic discontinuous Galerkin methods. InternationalJournal for Numerical Methods in Engineering (5), 752–784 (2006)33. Shu, C.W.: Different formulations of the discontinuous Galerkin method for the viscous terms. Z.-C. Shu, M. Mu, W. Xue, J.Zou (Eds.), Advances in Scientific Computing pp. 144–155 (2001)34. Sv¨ard, M., Nordstr¨om, J.: On the order of accuracy for difference approximations of initial-boundary value problems. Journalof Computational Physics (1), 333–352 (2006)35. Sv¨ard, M., Nordstr¨om, J.: On the convergence rates of energy-stable finite-difference schemes. Journal of Computational Physics , 108819 (2019)36. Worku, Z.A., Zingg, D.W.: Simultaneous approximation terms and functional accuracy for diffusion problems discretized withmultidimensional summation-by-parts operators. arXiv preprint arXiv:2012.07812 (2020)37. Yan, J., Crean, J., Hicken, J.E.: Interior penalties for summation-by-parts discretizations of linear second-order differentialequations. Journal of Scientific Computing75