We consider an exponentially fitted discontinuous Galerkin method for advection dominated problems and propose a block solver for the resulting linear systems. In the case of strong advection the solver is robust with respect to the advection direction and the number of unknowns

Ayuso, B.

Lombardi, Ariel

Pietra, Paola

Universitat Autònoma de Barcelona. Centre de Recerca Matemàtica

Zikatanov, Ludmil

Ayuso de Dios, Blanca

Centre de Recerca Matemàtica

Diposit Digital de Documents de la UAB

CRMPreprintSeriesnumber1055A BLOCK SOLVER FOR THE EXPONENTIALLY FITTEDIIPG-0 METHODBLANCA AYUSO DE DIOS, ARIEL LOMBARDI, PAOLA PIETRA,AND LUDMIL ZIKATANOVAbstract. We consider an exponentially fitted discontinuous Galerkinmethod for advection dominated problems and propose a block solver for theresulting linear systems. In the case of strong advection the solver is robustwith respect to the advection direction and the number of unknowns.1. IntroductionLet Ω ⊂ IR2 be a polygon, f ∈ L2(Ω), g ∈ H1/2(∂Ω) and let  > 0 be constant.We consider the advection-diffusion problem(1) −div(∇u− βu) = f in Ω, u = g on ∂Ω,where β ∈ [W 1,∞(Ω)]2 derives from a potential β = ∇ψ. In applications tosemiconductor devices, u is the electron density, ψ the electrostatic potential andthe electric field |∇ψ| might be fairly large in some parts of Ω, so that (1) be-comes advection dominated. Its robust numerical approximation and the designof efficient solvers, are still a challenge. Exponential fitting Brezzi et al. [2005]and discontinuous Galerkin (DG) are two approaches that have been combined inLombardi and Pietra [2011] to develop exponentially fitted DG methods (in pri-mal and mixed formulation). In this note, we consider a variant of these schemes,based on the use of the Incomplete Interior Penalty IIPG-0 method and proposean efficient solver for the resulting linear systems.The change of variable ρ := e−ψ u in the problem (1) leads to(2) −∇ · (κ∇ρ) = f in Ω, ρ = χ on ∂Ω ,where κ := eψ and χ := e−ψ g. An IIPG-0 approximation to (2) gives riseto the EF-IIPG-0 scheme for (1). We propose a block solver that uses ideasfrom Ayuso de Dios and Zikatanov [2009] and reduce the solution to that ofan exponentially fitted Crouziex-Raviart (CR) discretization, which has muchless degrees of freedom. The associated (CR) matrix is further reduced to anapproximate block lower triangular form, which is efficiently solved by a blockGauss-Siedel algorithm.In our description we focus on the case β = ∇ψ piecewise constant; although weinclude some numerical results for a more general case (cf. Test 2). Extensions of1CRMPreprintSeriesnumber10552 B. AYUSO DE DIOS, A. LOMBARDI, P. PIETRA, AND L. ZIKATANOVthe method (allowing ψ to be discontinuous) and further analysis of the proposedsolvers are topics of current research.2. The Exponentially Fitted IIPG-0 methodLet Th be a shape-regular family of partitions of Ω into triangles T and leth = maxT∈Th hT with hT denoting the diameter of T for each T ∈ Th. We assumeTh does not contain hanging nodes. We denote by Eoh and E∂h the sets of all interiorand boundary edges, respectively, and we set Eh = Eoh ∪ E∂h .Let T+ and T− be two neighboring elements, and n+, n− be their outwardnormal unit vectors, respectively (n± = nT±). Let ζ± and τ± be the restrictionof ζ and τ to T±. We define the average and jump trace operators:2{ζ} = (ζ+ + ζ−), [[ ζ ]] = ζ+n+ + ζ−n− on E ∈ Eoh,2{τ} = (τ+ + τ−), [[ τ ]] = τ+ · n+ + τ− · n− on E ∈ Eoh,and on e ∈ E∂h we set [[ ζ ]] = ζn and {τ} = τ . We will also use the notation(u,w)Th =∑T∈Th∫Tuwdx 〈u,w〉Eh =∑e∈Eh∫euwds ∀u,w,∈ V DG ,where V DG is the discontinuous linear finite element space defined by:V DG ={u ∈ L2(Ω) : u|T ∈ P1(T ) ∀T ∈ Th},Here, P1(T ) is the space of linear polynomials on T . Similarly, P0(T ) and P0(e)are the spaces of constant polynomials on T and e, respectively. For each e ∈ Eh,let P0e : L2(e) 7→ P0(e) (resp. P0T : L2(T ) 7→ P0(T ), for each T ∈ Th) be theL2-orthogonal projections defined byP0e (u) :=1|e|∫eu, ∀u ∈ L2(e) , P0T (v) :=1|T |∫Tv, ∀ v ∈ L2(T ) .We denote by V CR the classical Crouziex-Raviart (CR) space:V CR={v ∈ L2(Ω) : v|T ∈ P1(T )∀T ∈ Th and P0e [[ v ]] = 0 ∀ e ∈ Eh}.Note that v = 0 at the midpoint me of each e ∈ E∂h . To represent the functionsin V DG we use the basis {ϕe,T}T∈Th,e∈Eh , defined by(3) ∀T ∈ Th ϕe,T (x) ∈ P1(T ) e ⊂ ∂T ϕe,T (me′) = δe,e′ ∀e′ ∈ Eh .In particular, any w ∈ P1(T ) can be written as w = ∑e⊂∂T w(me)ϕe,T .We first consider the IIPG-0 approximation to the solution of (2): Find ρ ∈V DG such that A(ρ, w) = (f, w)Th for all w ∈ V DG with(4) A(ρ, w) = (κ∗T∇ρ,∇w)Th − 〈{κ∗T∇ρ}, [[w ]]〉Eh + 〈Se{[[ ρ ]]},P0([[w ]])〉Eh .CRMPreprintSeriesnumber1055BLOCK SOLVER FOR EXPONENTIAL FITTING IIPG-0 3Here, Se is the penalty parameter and κ∗T ∈ P0(T ) the harmonic average approx-imation to κ = eψ/ both defined in Lombardi and Pietra [2011] by:(5) κ∗T :=1P0T (κ−1)=P0T (e−ψ ), Se := αeh−1e {κ∗T}e .Next, following Lombardi and Pietra [2011] we introduce the local operatorT : V DG 7→ V DG that approximates the change of variable introduced before (2):(6) Tw :=∑T∈Th(Tw)|T =∑T∈Th∑e⊂∂TP0e (e−ψ )w(me)ϕe,T ∀w ∈ V DG .By setting ρ := Tu in (4), we finally get the EF-IIPG-0 approximation to (1):Find uh ∈ V DG s.t. B(uh, w) := A(Tuh, w) = (f, w)Th ∀w ∈ V DG with(7) B(u,w)=(κ∗T∇Tu,∇w)Th− 〈{κ∗T∇Tu}, [[w ]]〉Eh+ 〈Se{[[Tu ]]},P0[[w ]]〉Eh .It is important to emphasize that the use of harmonic average to approximate κ =eψ/ as defined in (5) together with the definition of the local approximation ofthe change of variables prevents possible overflows in the computations when |∇ψ|is large and  is small. (See Lombardi and Pietra [2011] for further discussion).Also, these two ingredients are essential to ensure that the resulting methodhas an automatic upwind mechanism built-in that allows for an accurate approx-imation of the solution of (1) in the advection dominated regime. We will discussthis in more detail in Section 3.Prior to close this section, we define for each e ∈ Eh and T ∈ Th:ψm,e := minx∈eψ(x) ψm,T := minx∈Tψ(x); ψm,T ≤ ψm,e for e ⊂ ∂T .In the advection dominated regime  |β|h = |∇ψ|h(8) P0T (e−(ψ/)) ' 2e−ψm,T P0ei(e−ψ/) '  e−ψm,e .The first of the above scalings together with the definitions in (5) implies(9) κ∗T '1eψm,T , Se ' α2|e|−1e(ψm,T1+ψm,T2) e = ∂T1 ∩ ∂T2 .3. Algebraic System & PropertiesLet A and B be the operators associated to the bilinear forms A(·, ·) (4) andB(·, ·) (7), respectively. We denote by A and B their matrix representation in thebasis {ϕe,T}T∈Th,e∈Eh (3). In this basis, the operator T defined in (6) is representedas a diagonal matrix, D, and B = AD. Thus, the approximation to (2) and (1)amounts to solve the linear systems (of dimension 2ne−nb; with ne and nb beingthe cardinality of Eh and E∂h , respectively):(10) Aρ = F , and Du = ρ or Bu = F˜ ,CRMPreprintSeriesnumber10554 B. AYUSO DE DIOS, A. LOMBARDI, P. PIETRA, AND L. ZIKATANOVwhere ρ,u,F and F˜ are the vector representations of ρ, u and the right handsides of the approximate problems. From the definition (6) of T it is easy todeduce the scaling of the entries of the diagonal matrix D = (di,i)2ne−nbi=1 .D = (di,j)2ne−nbi,j=1 di,i = P0ei(e−ψ/) '  e−ψm,e , di,j ≡ 0 i 6= j .We now revise a result from Ayuso de Dios and Zikatanov [2009]:Proposition 1. Let Z ⊂ V DG be the space defined byZ = {z ∈ L2(Ω) : z|T ∈ P1(T ) ∀T ∈ Th and P0e{v} = 0 ∀ e ∈ Eoh} .Then, for any w ∈ V DG there exists a unique wcr ∈ V CR and a unique wz ∈ Zsuch that w = wcr + wz , that is: V DG = V CR ⊕ Z. Moreover, A(wcr, wz) = 0∀wcr ∈ V CR, and ∀wz ∈ Z.Proposition 1 provides a simple change of basis from {ϕe,T} to canonical basisin V CR and Z that results in the following algebraic structure for (10):(11) ρ =[ρzρcr], A =[Azz 0Avz Avv], B =[Bzz 0Bvz Bvv].Due to the assumed continuity of ψ, D is still diagonal in this basis. The algebraicstructure (11) suggests the following exact solver:Algorithm 1. The solution u = uz + ucr satisfying B(u,w) = (f, w)Th, for allw ∈ V DG is then obtained by1. Solve for uz: B(uz, wz) = (f, wz)Th ∀wz ∈ Z.2. Solve for ucr: B(ucr, wcr) = (f, wcr)Th − B(uz, wcr) ∀wcr ∈ V CR.Next, wet discuss how to solve efficiently each of the above steps.Step 1: Solution in the Z-space. In Ayuso de Dios and Zikatanov [2009] itwas shown that Azz is a diagonal positive definite matrix. This is also true forBzz since it is the product of two diagonal matrices. The continuity of ψ implies(12) B(uz, wz) = 〈SeT[[uz ]],P0e ([[wz ]])〉Eh ∀ uz, wz ∈ Z .Using (8) and (5) we observe that the entries of Bzz scale as:Bzz = (bi,j)nei=1 bi,j = Sei |ei|djδi,j ' δi,jα2e−(ψm,e−ψm,T1−ψm,T2 )/which are always positive, so in particular Bzz it is also an M -matrix.Step 2: Solution in V CR. In Ayuso de Dios and Zikatanov [2009] it was shownthat the block Avv coincides with the stiffness matrix of a CR discretization of(2), and so it is an s.p.d. matrix. However, this is no longer true for Bvv whichis positive definite but non-symmetric.B(ucr, wcr) = (κ∗T∇Tucr,∇wcr)Th ∀ ucr , wcr ∈ V CR .CRMPreprintSeriesnumber1055BLOCK SOLVER FOR EXPONENTIAL FITTING IIPG-0 5In principle, the sparsity pattern of Bvv is that of a symmetric matrix. Using (8)and (5), we find that the entries of the matrix scale as:(13) Bvv =(bcri,j)ncr:=ne−nbi,jbcri,j := κ∗T|ei||ej||T | nei · nejdj ' e− (ψm,e−ψm,T ) .Since ψ is assumed to be piecewise linear, for each T , it attains its minimum (andalso its maximum) at a vertex of T , say x0 and ψm,e is attained at one of thevertex of the edge e, say xe. In particular, this implies thatψm,e − ψm,T ≈ ∇ψ · (xe − x0) = β · (xe − x0) ={0 xe = x0|β|h xe 6= x0Hence, in the advection dominated case  |β|h some of the entries in (13) vanish(up to machine precision) for  small; this is the automatic upwind mechanismintrinsic of the method. As a consequence, the sparsity pattern of Bvv is no longersymmetric and this can be exploited to re-order the unknowns so that Bvv canbe reduced to block lower triangular form.Notice also that for Th acute, the block Avv being the stiffness matrix of theCrouziex-Raviart approximation to (2), is an M-matrix. Hence, since the blockBvv is the product of a positive diagonal matrix and Avv, it will also be an M -matrix if the triangulation is acute (see Brezzi et al. [2005]).4. Block Gauss-Siedel solver for V CR-blockWe now consider re-orderings of the unknowns (dofs), which reduce Bvv toblock lower triangular form. For such reduction, we use the algorithm from Tarjan[1972] which roughly amounts to partitioning the set of dofs into non-overlappingblocks. In the strongly advection dominated case the size of the resulting blocksis small and a block Gauss-Seidel method is an efficient solver. Such techniqueshave been studied in Wang and Xu [1999] for conforming methods. The idea isto consider the directed graph G = (V ,E) associated with Bvv ∈ IRncr×ncr ; Ghas ncr vertices labeled V = {1, . . . , ncr} and its set of edges E has cardinalityequal to the number of nonzero entries1 of Bvv. By definition, (i, j) ∈ E iff bcrij 6=0. Note that in the advection dominated case, the built-in upwind mechanismresults in a non-symmetric sparsity pattern for Bvv (see the last two paragraphsof Section 3). Thus, we may have (i, j) ∈ E, while (j, i) /∈ E. Then, the problemof reducing Bvv to block lower triangular form of Bvv is equivalent to partitioningG as a union of strongly connected components. Such partitioning induces non-overlapping partitioning of the set of dofs, V = ∪Nbi=1ωi. For i = 1, . . . , Nb, letmi denote the cardinality of ωi; let Ii ∈ IRncr×mi be the matrix that is identityon dofs in ωi and zero otherwise; and Bvvi = ITi BvvIi is the block correspondingto the dofs in ωi. The block Gauss-Seidel algorithm reads: Let ucr0 be given, and1Each dof corresponds to a vertex in the graph; each nonzero entry to an edge.CRMPreprintSeriesnumber10556 B. AYUSO DE DIOS, A. LOMBARDI, P. PIETRA, AND L. ZIKATANOVassume ucrk has been obtained. Then ucrk+1 is computed via: For i = 1, . . . Nb(14) ucrk+i/Nb = ucrk+(i−1)/Nb + Ii(Bvvi )−1ITi(F − Bvvucrk+(i−1)/Nb).As we report in Section 5, the action of (Bvvi )−1 can be computed exactly sincein the advection dominated regime the size of the blocks Bvvi is small.5. Numerical ResultsWe present a set of numerical experiments to assess the performance of theproposed block solver. The tests refer to problem (2) with  = 10−3, 10−5, 10−7,and Ω is triangulated with a family of unstructured triangulations Th. In thetables given below J = 1 corresponds to the coarsest grid and each refined trian-gulation on level J , J = 2, 3, 4 is obtained by subdividing each of the T ∈ Th onlevel (J − 1) into four congruent triangles. From the number of triangles nT thetotal number of dofs for the DG approximation is 3nT .Test 1. Boundary Layer: Ω = (−1, 1)2, β = [1, 1]t, nT = 112 for the coarsestmesh and f is such that the exact solution is given byu(x, y) =(x+1 + e−2/ − 2e(x−1)/1− e−2/)(y +1 + e−2/ − 2e(y−1)/1− e−2/).Test 2. Rotating Flow: Ω = (−1, 1)×(0, 1), f = 0 and curlβ 6= 0,β =[2y(1− x2)−2x(1− y2)]tg(x, y) ={1 + tanh (10(2x+ 1)) x ≤ 0, y = 0,0 elsewhere .We stress that this test does not fit in the simple description given here, andspecial care is required (see Lombardi and Pietra [2011]). For the approximation,for each T ∈ Th, with barycenter (xT , yT ), we use the approximation β|T ≈∇ψ|T with ψ|T = 2yT (1− x2T )x− 2xT (1− 2y2T )y (and so ψ discontinuous). Thecoarsest grid has nT = 224 triangles. In Fig. 1 are represented the plot ofï1 ï0.8 ï0.6 ï0.4 ï0.2 0 0.2 0.4 0.6 0.8 1ï1ï0.8ï0.6ï0.4ï0.200.20.40.60.81Test 1:  50 blocks   156 Dofsï1 ï0.8 ï0.6 ï0.4 ï0.2 0 0.2 0.4 0.6 0.8 100.10.20.30.40.50.60.70.80.91Test 2:  122 blocks   316 DofsFigure 1. Plot of the connected components (blocks) of Bvv createdduring Tarjan’s algorithm: Test 1 with  = 10−5 (left); Test 2 with = 10−7 (right)CRMPreprintSeriesnumber1055BLOCK SOLVER FOR EXPONENTIAL FITTING IIPG-0 7the strongly connected components of the graph depicting the blocks for Bvvcreated during Tarjan’s algorithm, on the coarsest meshes; for Test 1 with  =10−5 (left figure) and for Test 2 with  = 10−7 (right figure). We have useddifferent line types (and colors) to distinguish strongly connected componentsin the directed graph. In Table 1 we report the number of blocks Nb createdduring Tarjan’s algorithm; the maximum size of the largest such block (Mb); theaverage block size (nav); and the number of block-Gauss-Seidel iterations. AfterTarjan’s algorithm is used to re-order the matrix Bvv, we use the block Gauss-Seidel algorithm (14) where each small block is solved exactly. In the tests thatwe report here and also in all other similar tests that we have done (with similaradvection dominance) the number of block-Gauss-Seidel iterations and the sizeof the blocks is uniformly bounded with respect to the number of dofs whenthe advection strongly dominates. Thus, the computational cost for one blockGauss-Seidel iteration in the advection dominated regime is the same as the costof performing a fixed number of matrix vector multiplications and the algorithmis optimal in such regime.AcknowledgmentsFirst author was partially supported by MEC grant MTM2008-03541. Secondauthor was supported by CONICET, Argentina. Last author was supported inpart by National Science Foundation NSF-DMS 0810982.ReferencesBlanca Ayuso de Dios and Ludmil Zikatanov. Uniformly convergent iterativemethods for discontinuous Galerkin discretizations. J. Sci. Comput., 40(1-3):4–36, 2009.F. Brezzi, L. D. Marini, S. Micheletti, P. Pietra, R. Sacco, and S. Wang. Dis-cretization of semiconductor device problems. I. In Handbook of numericalanalysis. Vol. XIII, pages 317–441. North-Holland, Amsterdam, 2005.Ariel Lombardi and P. Pietra. Exponentially fitted discontinuous galerkinschemes for singularly perturbed problems. numerical methods for partial dif-ferential equations. 2011. (to appear) doi: 10.1002/num.20701.Robert Tarjan. Depth-first search and linear graph algorithms. SIAM J. Comput.,1(2):146–160, 1972.Feng Wang and Jinchao Xu. A crosswind block iterative method for convection-dominated problems. SIAM J. Sci. Comput., 21(2):620–645, 1999.CRMPreprintSeriesnumber10558 B. AYUSO DE DIOS, A. LOMBARDI, P. PIETRA, AND L. ZIKATANOVTest 1HHHHHJ1 2 3 410−3Nb 44 150 484 1182Mb 23 47 95 191nav 3.55 4.32 5.45 9.02iters 7 19 43 16610−5Nb 50 210 866 3474Mb 23 47 95 191nav 3.12 3.08 3.05 3.07iters 4 4 4 1410−7Nb 50 210 866 3522Mb 23 47 95 191nav 3.12 3.08 3.05 3.03iters 4 4 4 4Test 2HHHHHJ1 2 3 410−3Nb 31 1 1 1Mb 211 1304 5296 21344nav 10.19 1304 5296 21344iters 10 1 1 110−5Nb 122 468 1822 7106Mb 4 4 7 37nav 2.59 2.78 2.91 3.00iters 4 4 7 2410−7Nb 122 468 1832 7247Mb 4 4 4 6nav 2.59 2.78 2.89 2.95iters 4 4 4 4Table 1. Number of blocks (Nb) created during the Tarjan’s or-dering algorithm, size of largest block (Mb), average size of blocks(nav) and number of block-Gauss-Seidel iterations (iters) for Test1 (left) and Test 2 (right) .CRMPreprintSeriesnumber1055BLOCK SOLVER FOR EXPONENTIAL FITTING IIPG-0 9B. Ayuso de DiosCentre de Recerca Matema`ticaBarcelona, SpainE-mail address: bayuso@crm.catA. LombardiDepartamento de Matema´ticaUniversidad de Buenos Aires & CONICETArgentinaE-mail address: aldoc7@dm.uba.arP. PietraIMATI-CNRPavia, ItalyE-mail address: pietra@imati.cnr.itL. ZikatanovDepartment of MathematicsPenn State UniversityUSAE-mail address: ltz@math.psu.edu

A block solver for the exponentally fitted IIPG-0 method

https://ddd.uab.cat/pub/prepub/2011/hdl_2072_203577/Pr1055.pdf

A block solver for the exponentally fitted IIPG-0 method

Abstract

Similar works

Full text

Available Versions

Diposit Digital de Documents de la UAB