Multigrid Methods for NIPG K. Johannsen The Center for Subsurface Modeling (CSM) The Institute for Computational Engineering and Sciences (ICES) The University of Texas, Austin, TX 78712, USA

Abstract In this paper we consider the class of NIPG discretizations for the Laplace equation on conforming, nondegenerate, quasiuniform grids in two space dimensions. We investigate different multigrid methods using block variants of the classical Jacobi and GaussSeidel type smoothers. Two classes of smoothing iterations are considered. The class of element block iterations appears to be efficient only for a restricted range of discretization schemes. In contrast, a new class of node-centered overlapping block iterations introduced in this paper exhibits a robust behavior with respect to the discretization penalty. The schemes are analyzed for a special one-dimensional model. In the two-dimensional setting thorough numerical experiments are carried out showing a robust behavior of the proposed methods. AMS Subject Classification: 65F10, 65F50, 65N22, 65N30, 65N55. Key words: smoothing property, multigrid method, sparse matrices, discontinuous Galerkin methods, robustness 1 Introduction The application of discontinuous Galerkin methods to elliptic boundary value problems began in the early 1970’s [1, 2, 4, 13, 22, 23, 33, 34]. Whereas the first methods were symmetric, in the late 1990’s the development of nonsymmetric discretizations gained increasing attention [3, 7, 12, 24–32]. To date, many aspects of the discontinuous

2

K. Johannsen

Galerkin methods are well understood. Due to their interesting properties, fast solution schemes for the resulting discrete equations became a field of active research. Especially attractive are multigrid methods because of their potential for optimal complexity. Their application to conforming finite element discretizations is well established. In the case of symmetric discontinuous Galerkin methods (SIPG [33] and LDG [11]) multigrid schemes have been applied successfully [9, 10, 14, 20, 21]. Nonsymmetric variants have been investigated in [6, 16–19]. In [6], a V-cycle multigrid iteration employing a block ILU-smoother is used to solve Poisson-type problems with varying coefficients discretized by Baumann’s method [3]. In numerical experiments the method is shown to be reasonably fast. Fourier analysis is applied to various discontinuous Galerkin methods investigating spectra and smoothing factors of the considered smoothing block iterations (Jacobi, Gauss-Seidel and symmetric Gauss-Seidel), as well as two-grid convergence rates in [16–18]. The analysis is restricted to the structured case. In [19], W-cycle convergence is shown for the NIPG method with an appropriately chosen penalty parameter. In this paper we discuss two different classes of block iterative schemes used as smoothers in multigrid. One corresponds to blocks associated with finite elements; this class of schemes has been analyzed previously. The one to be introduced here corresponds to a blocking of the degrees of freedom associated with the nodes of the grid. It has similarities with the smoothing iterations analyzed in [16–18]. Robustness issues regarding the resulting multigrid schemes with respect to the polynomial basis, the grids and the penalty parameter of the discretization are investigated. It is found that the proposed methods are far superior to the methods using element blocks, and they exhibit many of the desired robustness properties. The remainder of this paper is organized as follows. In section 2, we analyze two Jacobi-type iterations with respect to a one-dimensional model problem. In section 3, we define the two-dimensional model problem and describe the NIPG discretization. The multigrid methods are described in section 4. In section 5, we investigate the behavior of the multigrid methods as applied to the model problem defined in section 3. Finally, in section 6, concluding remarks are given. 2 Analysis of a one-dimensional model problem In this section we investigate the discretization of the one-dimensional Laplace equation on an infinite, equidistant grid using Baumann’s

Multigrid Methods for NIPG

3

method [3] with P2 -elements. We shall study the behavior of blocktype Jacobi smoothers. Let h > 0, eh,j = [jh, (j + 1)h], j ∈ Z and the infinite grid be defined by Th = {eh,j |j ∈ Z}. We associate with Th the broken Sobolev space Vh,2 = {u ∈ L2 |u|e ∈ P2 , e ∈ Th }, where P2 denotes the space of polynomials of order k ≤ 2. On V h,2 × Vh,2 we define the bilinear form XZ X X u0 v 0 dx − hu0 i[v] + hv 0 i[u], (1) ah (u, v) = e∈Th

e

x∈hZ

x∈hZ

where u0 denotes the derivative of u, hui(x) = 12 (u(x+ ) + u(x− )), with u(x+ ) = limy→x,y>x u(y) and u(x− ) correspondingly, and [u] = u(x+ ) − u(x− ). For details, see [3]. To represent vectors and linear mappings, we use the basis functions arising from the monomial basis {1, x, x2 } of P2 on the reference element er = [0, 1] and the appropriate affine transformation. Vh,2 is isomorphic to V2 = R3 ⊗ RZ , and linear mappings may be represented by infinite block matrices. Let the elements of Th be ordered lexicographically and     0 0 −6 0 1 2 1  1  1 1 1 ; 0 6 3 , D= L = −U T = 6h 2h 6 9 8 0 0 0 (1) is represented by an infinite block-tridiagonal Toeplitz matrix, which we write using the usual stencil notation A = [L|D|U ] ∈ V2 × V2 . We investigate A by means of Fourier analysis. Let k ∈] − π, π] and ϕk ∈ RZ with (ϕk )j = eikhj ,

j ∈ Z,

(2)

where i denotes the complex unit. Let x ∈ R 3 ; it follows that A(x ⊗ ϕk ) = (Ak x) ⊗ ϕk ,

(3)

with Ak = Le−ikh + D + U eikh ∈ C3×3 . 2.1 The element-block Jacobi iteration We investigate the block Jacobi iteration, obtained by blocking the degrees of freedom associated with each element. Let θ ∈ R; the damped variant is defined by S1 = I − θN1 A,

N1 = [0|D −1 |0],

4

K. Johannsen

where 0 denotes the 3 × 3 zero matrix and I the identity in V 2 . Like A, S1 is an infinite block-tridiagonal Toeplitz matrix. Let x ∈ R 3 , k ∈]-π, π], and ϕk ∈ RZ defined by (2); it holds that S1 (x ⊗ ϕk ) = (S1,k x) ⊗ ϕk ,

(4)

with S1,k = (1−θ)I−θD −1 Le−ikh −θD −1 U eikh ∈ C3×3 . In the present context we define the smoothing property of an iterative scheme as follows. Definition 1 Let k ∈] − π, π] and Ak , Sk ∈ C3×3 . Sk is said to fullfil the smoothing property if ν sup kAk S1,k k2 ≤ η(ν)

(5)

k∈]−π,π]

for a function η(ν) → 0, ν → ∞. Lemma 1 Let θ ∈]0, 2[ and Ak , S1,k defined according to (3), (4). It holds that lim

ν sup kAk S1,k k2 ≥ 1/3.

ν→∞ k∈]−π,π]

Proof We consider S1,k (θ), where we write its dependence on the damping parameter explicitly. Elementary manipulations show that   0 0 0 ν (6) A0 S1,0 (1) =  0 0 0  , ∀ν ≥ 1. 0 0 1/3 Let θ ∈]0, 2[. It holds that

S1,k (θ) = (1 − θ)I + θS1,k (1), and, therefore, ν   X ν

µ (1 − θ)ν−µ θ µ A0 S1,0 (1) µ µ=0 ν   X ν (1 − θ)ν−µ θ µ A0 S1,0 (1) + (1 − θ)ν A0 (I − S1,0 (1)) = µ

ν A0 S1,0 (θ)=

µ=0

=A0 S1,0 (1) + (1 − θ)ν A0 (I − S1,0 (1)).

Hence, ν ν sup kAk S1,k k2 ≥ kA0 S1,0 k2 ≥ 1/3 − (1 − θ)ν kA0 (I − S1,0 (1))k2 ,

k∈]−π,π]

which implies the initial assertion.

PSfrag replacements

ν ν kAk S2,k k2

Multigrid Methods for NIPG

5

2/ν Num. res.

1

0.1

0.01

0.001

1e-04 1

10

100

1000

Fig. 1. Numerical evaluation of the left hand side of (5) for k = 0, k = 10−α/100 π, α ∈ [0, 250]N and ν ∈ [1, 1000]N . The upper bound from (8) is indicated by the solid line.

The lemma shows that the smoothing property according to Definition 1 does not hold for S1 . From (6) we can see that the x2 component of the error is not reduced sufficiently for small Fourier frequencies k. Note that this behavior is independent of the damping factor. Due to these results, it is highly questionable that analogous multigrid schemes for higher dimensions (using Q 2 -elements) are applicable to Baumman’s scheme. A confirmation of this conclusion will be given in section 5.

2.2 An overlapping block Jacobi iteration An overlapping block Jacobi iteration is realized by an additive subspace correction method, where each subspace corresponds to the degrees of freedom associated with a pair of neighboring elements. Let θ ∈ R; this method is defined by S2 = I − θN2 A,

N2 = [NL |ND1 + ND2 |NU ],

6

K. Johannsen

where ND1 , ND2 , NL , NU ∈ R3×3 are given by    −1 N D1 N U DU = . N L N D2 LD Matrix S2 is an infinite, block-pentadiagonal Toeplitz matrix. Let x ∈ R3 , k ∈]-π, π], and ϕk ∈ RZ defined by (2); it holds that S2 (x ⊗ ϕk ) = (S2,k x) ⊗ ϕk

(7)

with an appropriate S2,k ∈ C3×3 . We will now investigate the smoothing property for this iteration. Lemma 2 Let θ = 1/2 and Ak , S2,k defined according to (3), (7). It holds that ν sup kAk S2,k k2 ≤ 2/ν,

∀ν ≥ 1.

(8)

k∈−π,π

Proof The proof shall not be given analytically but by numerical evaluation. The left hand side of (5) has been evaluated for k = 0, k = 10−α/100 π, α ∈ [0, 250]N and ν ∈ [1, 1000]N . The results are depicted in Fig. 1. They clearly confirm (8), and no reasonable doubts on the validity of assertion (8) are valid. Lemma 2 proves the smoothing property for the one-dimensional Laplace problem. Due to the strong assumptions made here it is not clear to what extent these results carry over to the two-dimensional situation. 3 Model problems and discretization Let Ω =]0, 1[2 ⊂ R2 and f ∈ L2 (Ω). Let u ∈ H01 (Ω) be the weak solution of −∆u = f, u = 0,

x ∈ Ω, x ∈ ∂Ω.

(9a) (9b)

The solution is fully regular, i.e. u ∈ H 2 (Ω) ∩ H01 (Ω). Let T be a con¯ forming, nondegenerate, quasiuniform triangulation with ∪ e∈T e = Ω, where e is a triangle or a quadrilateral. The nondegeneracy condition requires that there exists a parameter ρ > 0 such that, if he = diam(e), then e contains a ball with radius ρh e . With T we associate the broken Sobolev space H1 (T) = {u ∈ L2 (Ω)|u|e ∈ H1 (e)∀e ∈ T}.

Multigrid Methods for NIPG

7

Let ne denote the outer normal of e ∈ T and Γe,f := e¯∩ f¯ the common interface of e, f ∈ T. We define the average of a vector valued function u ∈ (H1 (T))d by ( u|e on ∂e ∩ ∂Ω hui = u|e + u|f )/2 on Γe,f ⊂ Γ \∂Ω and the jump of a scalar function u ∈ H 1 (T) by ( u|e ne on ∂e ∩ ∂Ω . [u] = u|e ne + u|f nf on Γe,f ⊂ Γ \∂Ω The variational formulation of (9) on grid T for the NIPG method reads as follows: Find u ∈ H1 (T) such that B(u, v)=L(v)

∀v ∈ H1 (T)

(10)

with B(u, v)=

XZ e∈T

∇u · ∇vdx −

Z

h∇ui · [v]ds

Γ

e

µ + h∇vi · [u]ds + h Γ Z

Z

[u] · [v]ds

(11)

Γ

and L(v)

=

XZ e∈T

f vdx. e

The parameter µ ∈ R is called the penalty parameter. We discretize the variational problem (10) by replacing H 1 (T) with the finite dimensional subspace VPk (T) = {u ∈ L2 (Ω)|u|e ∈ Pk (e)∀e ∈ T} ⊂ H1 (T), where Pk (e) denotes the space of complete polynomials on e with degree ≤ k. Since we are considering only element block type iterative schemes, the choice of the basis functions on each element is not relevant for the iterative schemes. Nevertheless, numerical stability issues and round-off error control suggest the usage of certain polynomials. We use L2 -orthogonalized basis functions. Alternatively, we use the space VQk (T) = {u ∈ L2 (Ω)|u|e ∈ Qk (e)∀e ∈ T} ⊂ H1 (T), where Qk (e) denotes the space of tensor product Ansatz functions with polynomials of degree up to k in each spatial direction.

8

K. Johannsen

4 Multigrid methods In this section we define the multigrid schemes that will be investigated subsequently. Let L ≥ 1, T0 be a conforming triangulation of Ω and Tl be the uniform refinement of Tl−1 , 0 < l ≤ L [5]. On each grid level we discretize the model problem (9) by means of the NIPG discretization. For the prolongation we use the natural embedding, for the restriction the L2 -projection. We use ν pre- and post-smoothing steps on each grid level and one recursive call for the coarse grid corrections (V (ν, ν)-cycle). The problem on the coarsest grid (level 0) is solved exactly. The components described so far are standard. For details we refer to [8, 15, 35, 37]. As smoothing iterations, we employ block variants of the Jacobi, the Gauss-Seidel and the symmetric Gauss-Seidel iterations. We shall describe them in the framework of subspace correction methods [36]. Let 0 < l ≤ L, k ≥ 1 and V be either VPk (Tl ) or VQk (Tl ). Let e ∈ Tl and Ve ⊂ V be the linear subspace spanned by the degrees of freedom of e. The solution u ∈ V of the discrete problem fulfills B(u, v)=L(v),

∀v ∈ V.

Let Pe : V → Ve be defined by B(Pe v, ve )=B(v, ve ),

∀v ∈ V, ve ∈ Ve

(12)

and u0 be an approximation of u. The subspace correction step corresponding to Ve leads to the new iterate un = u0 + Pe (u − u0 ). The correction ce = Pe (u − u0 ) corresponds to an exact solve on Ve , which can be found by solving B(ce , ve )=L(ve ) − B(u0 , ve ),

∀ve ∈ Ve

(13)

for ce ∈ Ve . Let θ ∈ R. The method of parallel subspace corrections X Pe (u − u0 ) (14) u1 = u 0 + θ e∈Tl

is referred to as the (damped) element block Jacobi iteration, [EBJAC]. It does not depend on the ordering of the elements. In order to define successive subspace correction methods, we must define an ordering of the subspaces. Let N = #T l and I : [1, N ]N → Tl be bijective. I induces an ordering of the elements and, subsequentially,

Multigrid Methods for NIPG Level Triangular grids Square grids

0 50 25

1 200 100

9 2 800 400

3 3.200 1.600

4 12.800 6.400

5 51.200 25.600

6 204.800 102.400

Table 1. Number of elements of the triangular grids (structured and unstructured) and of the square grids on different grid levels.

an ordering of the subspaces Ve , e ∈ T. We define the successive subspace correction method by u00 =u0 ,

(15a)

ui0 =u0 + Pei (u − ui−1 0 ),

0 < i ≤ N,

u1 =uN 0 .

(15b) (15c)

We refer to this iteration as the element block Gauss-Seidel iteration, [EB-GS]. The corresponding symmetric variant [EB-SGS] is defined by u00

=u0 , i u0 =u0 + +i uN =u0 + 0 u1 =u2N 0 .

(16a) Pei (u −

ui−1 0 ),

PeN −i+1 (u −

u0N +i−1 ),

0 < i ≤ N,

(16b)

0 < i ≤ N,

(16c) (16d)

Let the set of nodes of T be defined by ¯ = corner(e), e ∈ T}. N (T) = {p ∈ Ω|p For n ∈ N we define Vn = ⊕e∈T,e∩n6=∅ Ve , and Pn : V → Vn analogously to (12). The method of parallel subspace corrections, defined analogously to (14), is called the (damped) overlapping block Jacobi iteration, [OB-JAC]. Provided an ordering of the nodes is given, the overlapping block Gauss-Seidel, [OB-GS] and the symmetric, overlapping block Gauss-Seidel, [OB-SGS] are defined analogously to (15) and (16). 5 Numerical experiments In this section, we investigate the behavior of the multigrid methods described in the previous section. Since analytical results are not available we give a detailed discussion of numerical experiments

10

K. Johannsen

Fig. 2. Coarse grid triangulations used for the numerical experiments. The structured square grid (left) is investigated in section 5.1 and 5.2, the structured (middle) and unstructured (right) triangular grids are investigated in section 5.2.

in order to show the capabilities and the limits of the proposed schemes. For practical reasons, we limit our considerations to the two-dimensional model problem (9) on the unit square and to uniformly refined grid hierarchies. For the numerical experiments we use one of the grids shown in Fig. 2 as a coarse grid. The number of elements in the resulting grid hierarchies are given in Table 1. 5.1 Overlapping and non-overlapping iterations Now we compare the multigrid methods introduced in the previous section. The methods differ only in the smoothers employed. We investigate the schemes on the structured square grid hierarchy (arising from the coarse grid of Fig. 2, left) using the polynomial basis Q 2 on each element. This setting corresponds to a finite, tensor-product problem and resembles the situation analyzed in section 2, where P2 (= Q2 ) one-dimensional elements have been used. We investigate six V (ν, ν)-multigrid iterations using damped smoothers (damping factor θ) specified by: EB-JAC smoother, EB-GS smoother, EB-SGS smoother, OB-JAC smoother, OB-GS smoother, OB-SGS smoother,

θ = 0.75, θ = 1, θ = 1, θ = 0.5, θ = 1, θ = 1,

ν ν ν ν ν ν

= 4, = 1, = 1, = 4, = 1, = 1.

(17) (18) (19) (20) (21) (22)

The choice of the nontrivial damping factor for the Jacobi variants is experimental; the number of smoothing steps has been chosen to obtain reasonable convergence rates. (The ordering of the elements and nodes is chosen lexicographically; first left to right, then bottom to

(18) (19) (20) (21) (22) 1

1

10

0.01

Scheme (17) 100

1000

10000

100000

1e+06

1

0.001 0.001

PSfrag replacements Scheme (17) Scheme (18) Scheme (19) Scheme (20)

0.1

0.1

Scheme (22)

(19) (20) (21) (22) 0.01

PSfrag replacements Scheme (17)

Scheme Scheme Scheme Scheme 0.01

0.01

0.1

1

10

Scheme (20) 100

1000

0.1

1

10

PSfrag replacements Scheme (17) Scheme (18) Scheme (19) Scheme (20) Scheme (21)

0.01

PSfrag replacements Scheme (17) Scheme (18)

Scheme (20) Scheme (21) Scheme (22) 0.01

100000

1e+06

0.1

0.01

Scheme (18) 100

1000

10000

100000

1e+06

0.001 0.001

1

1

0.1

0.1

0.01

0.01

Scheme (21)

0.01

0.1

1

10

100

1000

0.01

0.1

1

10

100

1000

Scheme (19) 0.001 0.001

10000

1

0.1

0.001 0.001

11

1

0.1

0.001 0.001

PSfrag replacements Scheme (17) Scheme (18) Scheme (19)

Scheme (21) Scheme (22)

PSfrag replacements

Scheme Scheme Scheme Scheme Scheme

Multigrid Methods for NIPG

0.01

0.1

1

10

100

1000

10000

100000

10000

100000

1e+06

Scheme (22) 1e+06

0.001 0.001

10000

100000

1e+06

Fig. 3. Average convergence rates for the multigrid schemes (17)-(22). The calculations have been carried out on grid level l, l = 1, . . . , 5, corresponding to the different lines in the figures.

top. The geometrical position associated with the nodes is naturally its location; the position of the elements is its center of mass). We measure the convergence rate by   kd20 k2 1/10 ρ10,20 := , (23) kd10 k2 where di is the defect of the i-th iterate in a L 2 -orthogonal basis, defined by the right hand side of (13). We investigate the multigrid convergence for the different grid levels and penalty parameters l ∈ {1, . . . , 5},

µ ∈ {10−3 , 10−2 , . . . , 106 }.

(24)

The numerical results are given in Fig. 3. Each of the sub-figures corresponds to one of the multigrid schemes (17)-(22), as indicated, and displays the multigrid convergence rate in dependence on the

PSfrag replacements Scheme (17)? Scheme (18)? Scheme (19)?

(18)? (19)? (20)? (21)? (22)?

Scheme (21)? Scheme (22)?

PSfrag replacements

Scheme Scheme Scheme Scheme Scheme

12

1

K. Johannsen

1

Scheme (20)?

0.01

0.1

1

10

0.01

Scheme (17)? 100

1000

10000

100000

1e+06

1

0.001 0.001

PSfrag replacements Scheme (17)? Scheme (18)? Scheme (19)? Scheme (20)?

(19)? (20)? (21)? (22)? 0.001 0.001

0.1

Scheme (22)?

Scheme Scheme Scheme Scheme 0.01

PSfrag replacements Scheme (17)?

0.1

0.01

0.1

1

10

100

1000

10000

100000

1e+06

10

100

1000

10000

100000

1e+06

10

100

1000

10000

100000

1e+06

1

Scheme (21)?

0.001 0.001

0.01

0.1

1

10

PSfrag replacements Scheme (17)? Scheme (18)? Scheme (19)? Scheme (20)? Scheme (21)?

Scheme (20)? Scheme (21)? Scheme (22)? 0.01

PSfrag replacements Scheme (17)? Scheme (18)?

0.1

0.1

0.01

Scheme (18)? 100

1000

10000

100000

1e+06

1

0.001 0.001

0.01

0.1

1

1

Scheme (22)? 0.1

0.1

0.01

0.01

Scheme (19)? 0.001 0.001

0.01

0.1

1

10

100

1000

10000

100000

1e+06

0.001 0.001

0.01

0.1

1

Fig. 4. Average convergence rates for the multigrid schemes (17)-(22), where the number of smoothing steps is increased by a factor of four for each scheme (indicated by an asterisk). The calculations have been carried out on grid level l, l = 1, . . . , 5, corresponding to the different lines in the figures.

penalty parameter. The different graphs correspond to different grid levels (the graphs are not labeled, larger convergence rates correspond always to higher grid levels). In the left column the results corresponding to the element block iterations (from top to bottom: Jacobi, Gauss-Seidel, symmetric Gauss-Seidel) are given. They show mesh-size-independent convergence for a narrow band of penalty parameters, centered around 1. This result resembles the findings for the symmetric smoother investigated in [19] as well as the results obtained for the symmetric SIPG discretization reported in [14, 20]. As can be seen from the figure, the schemes do not converge for very small or very large penalty parameters. The results of the corresponding three overlapping block iterations are given in the right column of Fig. 3. In these cases the schemes exhibit mesh-size-independent

Multigrid Methods for NIPG

13

convergence even for small penalty parameters. Numerical results not given here indicate that this holds even for µ = 0, where the NIPG discretization coincides with Baumann’s method [3]. The behavior of the EB-JAC confirms, by analogy, the result for the one-dimensional case given in Lemma 1. The behavior of the OB-JAC confirms the result given in Lemma 2. For large penalties the schemes (17)-(22) exhibit increasingly poor multigrid convergence. For (17)-(19) this behavior cannot be improved, whereas for (20)-(22) it can be improved by increasing the number of smoothing steps. To illustrate this fact we display the convergence rates for (17)-(22) in Fig. 4, where we used a V (4ν, 4ν)-cycle instead of the V (ν, ν)-cycle investigated in Fig. 3. Whereas the modified schemes (17) ? -(19)? do not show a significant improvement in their large penalty behavior, the modified schemes (20)? -(22)? do (a further improvement of the modified scheme (20) may be obtained by further increasing the number of smoothing steps). The results in Fig. 3 and 4 show clearly that the usage of the overlapping block-type smoothers leads to multigrid schemes that are robust with respect to the penalty parameter. The results for the OB-JAC show that the robustness may be obtained by local smoothing procedures. Note that this method requires the choice of an appropriate damping factor.

5.2 The symmetric, overlapping block Gauss-Seidel iteration In this section we investigate the extent to which the results from the previous section carry over to more general cases. Previously, we investigated the behavior of different smoothers with respect to the model problem discretized by Q2 elements on regular square grids. It turned out that the overlapping block-type smoothers show robustness with respect to the penalty parameter. Since we are interested in the development of practically applicable multigrid schemes we study the scheme using the OB-SGS (22). We investigate its convergence properties for different grid hierarchies arising form the coarse grids given in Fig. 2 and for the different element Ansatz spaces P k and Qk , k = 1, . . . , 5. We consider the same range of penalty parameters as in the previous section (24). The range of refinement levels is restricted by memory limitations, and it depends on the element Ansatz spaces. The maximum levels investigated are lP1 = 6, lQ1 = 6,

lP2 = 6, lQ2 = 5,

lP3 = 5, lQ3 = 4,

lP4 = 4, lQ4 = 4,

lP5 = 4, lQ5 = 3.

14

K. Johannsen

5.2.1 Structured square grids The results for the square grids are given in Fig. 5. In the left column the convergence rates for P k , k = 1, . . . , 5 are displayed and in the right coulumn are rates for Q k , k = 1, . . . , 5. In the discussion of our results, we distinguish three regions of the penalty parameter. Small penalties (0 ≤ µ  1) This range of penalties resembles the behavior of its limit at µ = 0, Baumann’s scheme. It can be seen clearly from Fig. 5 that the schemes using P 1 or Q1 elements lead to the divergence of the multigrid method. This is due to a lack of approximation in the resulting discretization. For the P k -, Qk -elements, where k ≥ 2, an L2 -convergence of second order is obtained [27], and the approximation is sufficient for mesh-size-independent multigrid convergence. In Fig. 5 h-uniform convergence is documented, indicating that the OB-SGS iteration is a smoother for any of the discretizations considered. By analogy, the smoothing properties are explained by the findings given in section 2. Medium penalties (µ ≈ 1) In this range of penalties, all discretization schemes have good approximation qualities (second order in L 2 ). For polynomial degree k = 1 there is numerical evidence; for k ≥ 2 a proof has been given [27]. The OB-SGS shows good smoothing qualities and an h-independent multigrid convergence result in all cases. By analogy, the smoothing properties observed here may be explained by the behavior seen in the case of the symmetric SIPG discretization. In that case, even the EB-SGS has been shown to obey the smoothing property [14]. For the nonsymmetric case considered here, a variant of the Gauss-Seidel iteration has been shown to be a smoother [19]. Large penalties (µ  1) In this case we shall distinguish different situations according to the properties of the penalty term from (11), i.e. Z 1 [u] · [v]ds. (25) h Γ For the discretization scheme to exhibit an sufficient L 2 -convergence for increasing penalties, the kernel of (25) must be sufficiently large. For P1 elements it is easy to see that the kernel is trivial. Mesh size independent multigrid convergence therefore cannot be expected. This is in accordance with the numerical findings. In all other cases the kernel contains at least the conforming Q 1 -functions and is therefore large enough for a sufficient approximation. For P k , k ≥ 2, the

Multigrid Methods for NIPG

15

numerical experiments show a deterioration of the multigrid convergence with increasing penalty, indicating a failure of the smoothing iteration. In the case of tensor product elements (Q k , k ≥ 1), the numerical results for the V (1, 1)-cycle show a slight increase in the multigrid convergence rate with increasing penalty. Numerical experiments not documented here show that the usage of a W (ν, ν)cycle multigrid leads to h-independent multigrid convergence for Q k elements, provided ν is large enough. This does not hold true in the case of Pk -elements. An explanation of this behavior is beyond the scope of this paper. 5.2.2 Structured triangular grids The results for the structured triangular grids are given in Fig. 6. The behavior of the multigrid convergence for small and medium penalty parameters resembles the one found for the square grids. The argumentation is the same, even though the similarity with the investigations made in section 2 is not as striking as in the previous case. The major differences can be seen for large penalties. For Pk -elements, k ≥ 1, the kernel of (25) contains the conforming P1 -functions, i.e. the kernel is large enough to obtain a sufficient approximation. This is in accordance with the mesh-size-independent multigrid convergence observed for P k , k ≥ 1. The numerical results for Qk , k ≥ 1 show a deterioration of the multigrid convergence with increasing penalty parameter. Results not documented here show that for k = 1, 2, 3 the usage of a W (ν, ν)-cycle, with ν large enough, is sufficient for fast multigrid convergence in the large penalty range. For k = 4, 5 a remedy of this type is not possible. An explanation of this behavior is beyond the scope of this paper. 5.2.3 Unstructured triangular grids The results for the unstructured grid are given in Fig. 7. They resemble the behavior for structured triangular grids, indicating that the structured nature of the grid considered in the previous section is of no importance to the multigrid convergence behavior. 6 Conclusions In this paper we have investigated several variants of the classical iterative schemes (Jacobi, Gauss-Seidel, symmetric Gauss-Seidel) as applied to the NIPG discontinuous Galerkin discretization of the Laplace equation. The variants may be divided into two groups. On the one hand we considered the methods arising from a blocking of the unknown associated with the elements. These methods lead to

16

K. Johannsen

multigrid methods with h-independent convergence only for penalty parameters within a narrow band. Robust behavior with respect to the penalization cannot be obtained. The methods are especially not applicable to Baumann’s methods with vanishing penalty. On the other hand we investigated the behavior of overlapping block schemes, where the unknowns are grouped by association to the nodes of the (conforming) underlying grid. These methods are better described in terms of subspace correction methods. In many cases, they lead to robust behavior with respect to the penalty. They are especially well-suited for Baumann’s methods with vanishing penalty. A first theoretical approach has been presented, explaining the behavior of these two groups of schemes. Acknowledgments The support of this work by the J. Tinsley Oden Faculty Fellowship Research Fund is gratefully acknowledged. Furthermore, the support of P. Bastian, Heidelberg University, Germany, who made a software implementation of discontinuous Galerkin methods available is gratefully acknowledged. References 1. D.N. Arnold. An iterior penalty finite element method with discontinuous elements. J. Numer. Anal., 19:742–760, 1982. 2. I. Babuska and M. Zlamal. Nonconforming elements in the finite element method with penalty. SIAM J Numer. Anal., 10:863–875, 1973. 3. I. Babusky, C.E. Baumann, and J.T. Oden. A discontinuous hp finite element method for diffusion problems. J. Comp. Phys., 146:491–516, 1998. 4. G. Baker. Finite elment methods for elliptic equations using nonconforming elements. Math. Comp., 31:45–59, 1977. 5. R.E. Bank, A.H. Sherman, and A. Weiser. Refinement algorithms and data structures for regular local mesh refinement. In: Scientific computing IMACS, North-Holland, Amsterdam, 1983. 6. P. Bastian and V. Reichenberger. Multigrid for higher order discontinuous Galerkin finite elements applied to groundwater flow. Technical Report 200037, SFB 359, Heidelberg University, 2000. 7. C.E. Baumann and J.T. Oden. A discontinuous hp finite element method for convection-diffusion problems. Comp. Meth. Appl. Mech. Eng., 175(3-4):311– 341, 1999. 8. J.H. Bramble and J.E. Pasciak. New convergence estimates for multigrid algorithms. Math. Comp., 49:311–329, 1987. 9. S. Brenner and Li-Yeng Sung. c0 interior penalty methods for fourth order elliptic boundary value problems on polygonal domains. J. Sci. Comp., 2223:83–118, 2005.

Multigrid Methods for NIPG

17

10. S. Brenner and Jie Zhao. Convergence of multigrid algorithms for interior penalty methods. Appl. Num. Anal. Comp. Math., 2(1):3–18, 2005. 11. B. Cockburn and C.W. Shu. The local discontinuous Galerkin finite element method for convection diffusion systems. SIAM J. Num. Anal., 35:2440–2463, 1998. 12. C. Dawson, S. Sun, and M.F. Wheeler. Compatible algorithms for coupled flow and transport. Comp. Meth. Appl. Mech. Eng., 193:2565–2580, 2004. 13. J. Douglas and T. Dupont. Interior penalty procedures for elliptic and parabolic Galerkin methods. Lecture notes in Physics, 58:207–216, 1976. 14. J. Gopalakrishnan and G. Kanschat. A multilevel discontinuous Galerkin method. Numer. Math., 3:527–550, 2003. 15. W. Hackbusch. Multi-Grid Methods and Applications. Springer-Verlag, Berlin, Heidelberg, 1985. 16. P.W. Hemker, W. Hoffmann, and M.H. van Raalte. Fourier two-level analysis for discontinuous Galerkin discretization with linear elements. Technical Report MAS-R0217, ISSN 1386-3703, CWI, Amsterdam, 2002. 17. P.W. Hemker, W. Hoffmann, and M.H. van Raalte. Two-level fourier analysis of a multigrid approach for discontinuous Galerkin discretisation. Technical Report MAS-R0206, ISSN 1386-3703, CWI, Amsterdam, 2002. 18. P.W. Hemker and M.H. van Raalte. Fourier two-level analysis for higher dimensional discontinuous Galerkin. Technical Report MAS-R0227, ISSN 1386-3703, CWI, Amsterdam, 2002. 19. K. Johannsen. A symmetric smoother for the nonsymmetric interior penalty discontinuous galerkin discretization. Technical Report ICES Report 05-23, University of Texas at Austin, USA, 2005. 20. G. Kanschat. Preconditioning methods for local discontinuous Galerkin discretizations. SIAM J. Sci. Comp., 25(3):815–831, 2003. 21. G. Kanschat. Block preconditioner for ldg discretizations of linear incompressible flow problems. J. Sci. Comp., 22-23:371–384, 2005. 22. J.T. Oden and L.C. Wellford Jr. Discontinuous finite element approximations for the analysis of shock waves in nonlinear elastic materials. J. Comp. Phys., 19(2):179–210, 1975. 23. P. Percell and M.F. Wheeler. A local residual finite element procedure for elliptic equations. SIAM J. Num. Anal., 15(4):705–714, 1978. 24. B. Riviere. Discontinous Galerkin finite element methods for solving the miscible displacement problem in porous media. PhD thesis, The University of Texas at Austin, 2000. 25. B. Riviere and M.F. Wheeler. Non conforming methods for transport with nonlinear reaction. Contemporary Mathematics, 295:421–432, 2002. 26. B. Riviere, M.F. Wheeler, and V. Girault. Improved energy estimates for inerior penalty, constraint and discontinuous Galerkin methods for elliptic problems. Part i. Computational Geosciences, 3:337–360, 1999. 27. B. Riviere, M.F. Wheeler, and V. Girault. A priori error estimates for finite element methods based on discontinuous approximation spaces for elliptic problems. SIAM J. Numer. Anal., 39(3):902–931, 2001. 28. S. Sun. Discontinuous Galerkin methods for reactive transport in porous media. PhD thesis, The University of Texas at Austin, 2003. 29. S. Sun and M.F. Wheeler. Discontinuous Galerkin methods for coupled flow and reactive transport problems. Appl. Num. Math., 52:273–298, 2005.

18

K. Johannsen

30. S. Sun and M.F. Wheeler. L2 (H 1 ) norm a posteriori error estimation for discontinuous Galerkin approximations of reactive transport problems. J. Sci. Comp., 22:501–530, 2005. 31. S. Sun and M.F. Wheeler. Anisotropic and dynamic mesh adaptation for discontinuous Galerkin methods applied to reactive transport. Comput. Meth. Appl. Mech. Eng., to appear. 32. S. Sun and M.F. Wheeler. Symmetric and non-symmetric discontinuous Galerkin methods for reactive transport in porous media. SIAM J. Num. Anal., to appear. 33. M.F. Wheeler. An elliptic collocation-finite element method with interior penalties. SIAM J Numer. Anal., 15:152–161, 1978. 34. M.F. Wheeler and B.L. Darlow. Interior penalty Galerkin procedures for miscible displacement problems in porous media. Comp. Meth. Nonlin. Mech., pages 485–506, 1980. 35. J. Xu. Iterative methods by space decomposition and subspace correction. SIAM Rev., 34(4):581–613, 1992. 36. J. Xu and L. Zikatanov. The method of alternating projections and the method of subspace corrections in hilbert space. J. Amer. Math. Soc., 15:573– 597, 2002. 37. H. Yserentant. Old and new convergence proofs for multigrid methods. Acta Numerica, pages 285 – 326, 1993.

PSfrag replacements p1 p2 p3 p4 p5

PSfrag replacements

Multigrid Methods for NIPG

q2 q3 q4 q5

p2 p3 p4 p5 q1 q2 q3 q4 q5 1

19

1

p1

q1

0.001 0.001

0.01

0.1

1

10

100

1000

10000

0.1

0.01

100000

1e+06

0.001 0.001

0.01

q3 q4 q5

p3 p4 p5 q1 q2 q3 q4 q5 1

PSfrag replacements p1 p2 p3 p4 p5 q1

0.01

PSfrag replacements p1

0.1

0.1

1

10

100

1000

10000

100000

1e+06

100000

1e+06

100000

1e+06

1

p2

q2

0.001 0.001

0.01

0.1

1

10

100

1000

10000

0.1

0.01

100000

1e+06

0.001 0.001

0.01

q4 q5

p4 p5 q1 q2 q3 q4 q5 1

PSfrag replacements p1 p2 p3 p4 p5 q1 q2

0.01

PSfrag replacements p1 p2

0.1

0.1

1

10

100

1000

10000

1

p3

q3

0.001 0.001

0.01

0.1

1

10

100

1000

10000

0.1

0.01

100000

1e+06

0.001 0.001

q5

p5 q1 q2 q3 q4 q5 1

PSfrag replacements p1 p2 p3 p4 p5 q1 q2 q3

0.01

PSfrag replacements p1 p2 p3

0.1

0.01

0.1

1

10

100

1000

10000

1

p4

q4

0.001 0.001

0.01

0.1

1

10

100

1000

10000

q1 q2 q3 q4 q5

PSfrag replacements p1 p2 p3 p4 p5 q1 q2 q3 q4

0.01

PSfrag replacements p1 p2 p3 p4

0.1

0.1

0.01

100000

1e+06

1

0.001 0.001

0.01

0.1

1

10

100

1000

10000

100000

1e+06

0.01

0.1

1

10

100

1000

10000

100000

1e+06

1

p5

q5

0.1

0.1

0.01

0.01

0.001 0.001

0.01

0.1

1

10

100

1000

10000

100000

1e+06

0.001 0.001

Fig. 5. Average convergence rates for V (1, 1)-cycle multigrid method using the symmetric, overlapping block Gauss-Seidel iteration as smoother on a structured mesh hierarchy using square elements. The left column corresponds to the p1 , . . . , p5 basis (top to bottom). The right column corresponds to the q1 , . . . , q5 basis (top to bottom).

K. Johannsen

q2 q3 q4 q5

p2 p3 p4 p5 q1 q2 q3 q4 q5

PSfrag replacements p1 p2 p3 p4 p5

PSfrag replacements

20 1

1

p1

q1

0.001 0.001

0.01

0.1

1

10

100

1000

10000

0.1

0.01

100000

1e+06

0.001 0.001

0.01

q3 q4 q5

p3 p4 p5 q1 q2 q3 q4 q5 1

PSfrag replacements p1 p2 p3 p4 p5 q1

0.01

PSfrag replacements p1

0.1

0.1

1

10

100

1000

10000

100000

1e+06

100000

1e+06

100000

1e+06

1

p2

q2

0.001 0.001

0.01

0.1

1

10

100

1000

10000

0.1

0.01

100000

1e+06

0.001 0.001

0.01

q4 q5

p4 p5 q1 q2 q3 q4 q5 1

PSfrag replacements p1 p2 p3 p4 p5 q1 q2

0.01

PSfrag replacements p1 p2

0.1

0.1

1

10

100

1000

10000

1

p3

q3

0.001 0.001

0.01

0.1

1

10

100

1000

10000

0.1

0.01

100000

1e+06

0.001 0.001

q5

p5 q1 q2 q3 q4 q5 1

PSfrag replacements p1 p2 p3 p4 p5 q1 q2 q3

0.01

PSfrag replacements p1 p2 p3

0.1

0.01

0.1

1

10

100

1000

10000

1

p4

q4

0.001 0.001

0.01

0.1

1

10

100

1000

10000

q1 q2 q3 q4 q5

PSfrag replacements p1 p2 p3 p4 p5 q1 q2 q3 q4

0.01

PSfrag replacements p1 p2 p3 p4

0.1

0.1

0.01

100000

1e+06

1

0.001 0.001

0.01

0.1

1

10

100

1000

10000

100000

1e+06

0.01

0.1

1

10

100

1000

10000

100000

1e+06

1

p5

q5

0.1

0.1

0.01

0.01

0.001 0.001

0.01

0.1

1

10

100

1000

10000

100000

1e+06

0.001 0.001

Fig. 6. Average convergence rates for V (1, 1)-cycle multigrid method using the symmetric, overlapping block Gauss-Seidel iteration as smoother on a structured mesh hierarchy using triangular elements. The left column corresponds to the p1 , . . . , p5 basis (top to bottom). The right column corresponds to the q1 , . . . , q5 basis (top to bottom).

PSfrag replacements p1 p2 p3 p4 p5

PSfrag replacements

Multigrid Methods for NIPG

q2 q3 q4 q5

p2 p3 p4 p5 q1 q2 q3 q4 q5 1

21

1

p1

q1

0.001 0.001

0.01

0.1

1

10

100

1000

10000

0.1

0.01

100000

1e+06

0.001 0.001

0.01

q3 q4 q5

p3 p4 p5 q1 q2 q3 q4 q5 1

PSfrag replacements p1 p2 p3 p4 p5 q1

0.01

PSfrag replacements p1

0.1

0.1

1

10

100

1000

10000

100000

1e+06

100000

1e+06

100000

1e+06

1

p2

q2

0.001 0.001

0.01

0.1

1

10

100

1000

10000

0.1

0.01

100000

1e+06

0.001 0.001

0.01

q4 q5

p4 p5 q1 q2 q3 q4 q5 1

PSfrag replacements p1 p2 p3 p4 p5 q1 q2

0.01

PSfrag replacements p1 p2

0.1

0.1

1

10

100

1000

10000

1

p3

q3

0.001 0.001

0.01

0.1

1

10

100

1000

10000

0.1

0.01

100000

1e+06

0.001 0.001

q5

p5 q1 q2 q3 q4 q5 1

PSfrag replacements p1 p2 p3 p4 p5 q1 q2 q3

0.01

PSfrag replacements p1 p2 p3

0.1

0.01

0.1

1

10

100

1000

10000

1

p4

q4

0.001 0.001

0.01

0.1

1

10

100

1000

10000

q1 q2 q3 q4 q5

PSfrag replacements p1 p2 p3 p4 p5 q1 q2 q3 q4

0.01

PSfrag replacements p1 p2 p3 p4

0.1

0.1

0.01

100000

1e+06

1

0.001 0.001

0.01

0.1

1

10

100

1000

10000

100000

1e+06

0.01

0.1

1

10

100

1000

10000

100000

1e+06

1

p5

q5

0.1

0.1

0.01

0.01

0.001 0.001

0.01

0.1

1

10

100

1000

10000

100000

1e+06

0.001 0.001

Fig. 7. Average convergence rates for V (1, 1)-cycle multigrid method using the symmetric, overlapping block Gauss-Seidel iteration as smoother on an unstructured mesh hierarchy using triangular elements. The left column corresponds to the p1 , . . . , p5 basis (top to bottom). The right column corresponds to the q1 , . . . , q5 basis (top to bottom).

Multigrid Methods for NIPG

The Center for Subsurface Modeling (CSM). The Institute for Computational Engineering and Sciences (ICES). The University of Texas, Austin, TX 78712, USA. Abstract In this paper we consider the class of NIPG discretiza- tions for the Laplace equation on conforming, nondegenerate, quasi- uniform grids in two space ...

601KB Sizes 0 Downloads 202 Views

Recommend Documents

Multigrid Methods with Space-Time Concurrency - Computation
Key words. multigrid methods; space-time discretizations; parallel-in-time ... computing leading towards computers with more, but not faster, processors induce ...... inverse bandwidth cost, i.e., the cost per amount of data sent, and γ is the flop 

Multigrid methods with space-time concurrency - Computation
resources than standard space-parallel methods with se- quential time stepping. ...... Friedhoff, S., MacLachlan, S.: A generalized predictive analysis tool for ...

Algebraic Multigrid Methods for Laplacians of Graphs
Then the Galerkin operator ˆA ∈ Rm×m is defined as. ˆA = PT AP. ...... behind these goals comes from an interpretation in parallel computing. The subgraphs.

Multigrid methods with space-time concurrency
Department of Computer Science, KU Leuven, .... of the degree of anisotropy in the discrete operator. ... Relaxation is accelerated using a coarse-grid correc-.

Multigrid methods with space-time concurrency
E-mail: [email protected] ... plications including fluid flow, magnetohydrodynamics, compressible flow, and charged particle transport. Cur-.

A multigrid-in-time algorithm for solving evolution ...
our multigrid-in-time algorithm simply calls an existing time-stepping routine. However, to ...... Algebraic Multigrid Cycle on HPC Platforms, in 25th ACM International Conference on Supercomputing,. Tucson, AZ ... App. Math. and Comp. Sci., 5.

A parallel multigrid Poisson solver for fluids simulation ...
We present a highly efficient numerical solver for the Poisson equation on irregular voxelized domains ... a preconditioner for the conjugate gradient method, which enables the use of a lightweight, purely geometric ..... for transferring data across

An experimental investigation of Algebraic Multigrid for ...
Sep 9, 2011 - 'broadband' correction, one that ideally will address all components of the error spectrum .... complexity configurations are usually much more cost effective; hence, the moderate complexity ...... 10 and 11 (compare with. Figs.

An Algebraic Multigrid Method for the Stokes Equations
densation leading to a diffusion type term for the pressure-pressure coupling. The degrees of freedom are therefore collocated in the nodes of a conforming grid.

Parallel time integration with multigrid
In the case that f is a linear function of u(t), the solution to (1) is defined via ... scalability for cases where the “coarse-in-time” grid is still too large to be treated ...

Automated Methods for Evolutionary Pavé Jewellery Design
Jan 15, 2006 - Keywords Jewellery design, evolutionary algorithm, aesthetics, ..... Whilst the more natural application of this algorithm might appear to be in ...... to aid in the automated construction of Pavé jewellery exist, although at a price.

NUMERICAL METHODS FOR UNCERTAINTY ...
A straightforward application of Monte Carlo (MC) method leads to the .... obviously not the case in the PC approach that requires the development of a modified.