Repositorio Institucional de Documentos

Trabajo fin de máster Máster en Mecánica aplicada Programa oficial de postgrado en Mecánica computacional Curso 2011-2012 Métodos multimalla geométricos en mallas semi-estructuradas de Voronoi Pablo Salinas Cortés 5 de junio de 2012 Director: Francisco José Gaspar Lorenz Departamento de matemática aplicada Escuela de Ingeniería y Arquitectura Universidad de Zaragoza Métodos multimalla geométricos en mallas semi-estructuradas de Voronoi: Resumen Las ecuaciones en derivadas parciales son ampliamente utilizadas para modelizar gran cantidad de problemas físicos, debido a que son el tipo de ecuaciones que se usa para representar la naturaleza. Sin embargo, aunque están ampliamente extendidas en su uso, no siempre se pueden resolver analíticamente, y hay por tanto, que resolverlas numéricamente discretizándolas en una malla. Esta metodología, permite obtener los resultados deseados dentro de un dominio dadas unas condiciones de contorno. Por contra, los sistemas de ecuaciones obtenidos son muy grandes. Esto conlleva escoger muy cuidadosamente qué metodología se va a utilizar para resolverlo si no queremos malgastar recursos y tiempo. Dentro de los posibles algoritmos, el método multimalla destaca por su velocidad y por que el número de iteraciones necesarias para resolver es independiente del número de incógnitas. Dentro de los métodos multimalla, nos encontramos dos, el geométrico y el algebraico. El primero, es muy rápido pero solo se puede aplicar en dominios simples, el algebraico por otra parte es más lento pero se puede usar en formas complejas. Nosotros, en este trabajo hemos utilizado un método semi-estructurado que consiste en realizar una triangulación no estructurada sobre el dominio, como en el caso algebraico, para captar bien la superficie, y posteriormente dentro de cada triángulo, subdividirlos geométricamente, conectando los puntos medios de los lados cuantas veces sea necesario para obtener la precisión deseada. Además, dentro de cada triángulo de la malla no estructurada se pueden seguir diferentes estrategias dependiendo de las características de este, por ejemplo si es equilátero o isósceles, permitiendo una mayor optimización. Concretamente, hemos utilizado una malla de Delaunay como triangulación grosera y posteriormente hemos discretizado mediante volúmenes finitos. Como centros de los triángulos hemos usado los centros de Voronoi, toda malla Delaunay tiene asociada una de Voronoi, por sus buenas características; como por ejemplo que este centro es siempre perpendicular a los lados, lo que permite aproximar las derivadas con diferencias finitas. Sin embargo, el uso de estos puntos, a pesar de sus ventajas, introduce anisotropía en la discretización, exige el uso de triángulos acutángulos, así como al ir avanzando a lo largo de la jerarquía de mallas los puntos no están anidados, es decir, al restringir o prolongar los puntos no coinciden en coordenadas. Esto significa, que el método multimalla va a perder rendimiento si no hacemos nada para compensarlo. Para solucionar estos dos problemas hemos diseñado un conjunto de suavizadores que compensan los problemas de anisotropía y de no anidamiento, obteniendo unos resultados muy satisfactorios. Fijándonos en casos más particulares, podemos tener materiales cuyas propiedades varíen de un punto a otro, por ejemplo, un material compuesto. Si la variación de estas propiedades es de varios órdenes de magnitud, el algoritmo multimalla convencional no va a funcionar correctamente. Sin embargo, este problema puede solventarse si en la malla más fina se utiliza la media harmónica en las zonas de contacto y en las mallas bastas del multimalla se tiene en cuenta que se debe mantener la continuidad en el flujo. Esto se consigue utilizando el operador de Galerkin en las mallas más bastas. A su vez, hemos propuesto un algoritmo que permite obtener este operador en el caso de mallas triangulares, utilizando solamente moléculas en su cálculo. Los métodos propuestos, tanto los suavizadores, como el operador de Galerkin han sido puestos a prueba en varios ejemplos numéricos para comprobar su eficacia. Contents 1 Introduction 2 2 Introduction to Multigrid 5 3 Discretization of a diffusion problem on Voronoi grids 9 3.1 Discretization on triangular unstructured grids . . . . . . . . . . 9 3.2 Discretization on triangular structured grids . . . . . . . . . . . . 10 3.2.1 Stencil dependence on two angles characterizing the triangulargrid. ......................... 13 4 Multigrid method 15 4.1 Coarse-grid correction . . . . . . . . . . . . . . . . . . . . . . . . 15 4.2 Smoothers............................... 17 4.2.1 Jacobi smoother. . . . . . . . . . . . . . . . . . . . . . . . 17 4.2.2 Red-Black smoother. . . . . . . . . . . . . . . . . . . . . . 17 4.2.3 Diamond smoothers. . . . . . . . . . . . . . . . . . . . . . 18 4.2.4 Wormy smoothers. . . . . . . . . . . . . . . . . . . . . . . 19 5 Results of the proposed multigrid method on structured grids 20 6 Numerical experiments 24 6.1 Laplace problem in an A-shaped domain . . . . . . . . . . . . . . 25 6.2 Convection-diffusion problem on a square domain . . . . . . . . . 25 7 Piecewise discontinuous diffusion coefficient 28 7.1 Galerkinoperator........................... 29 7.2 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . 31 7.2.1 Diffusion problem on the unit square with various diffusioncoefficients........................ 31 7.2.2 Diffusion problem on a composite material domain . . . . 33 8 Conclusions 35 1 Chapter 1 Introduction In this work, we are interested in the multigrid solution of the large sparse system of equations, arising from the cell-centered finite volume discretization of a twodimensional partial differential equation on acute Delaunay triangulations. The use of Delaunay grids is motivated by the good properties of the discretizations obtained on this type of triangulations; by requiring them to be acute, we achieve monotonicity properties, see [4, 5, 16]. Every Delaunay grid has an associated dual mesh, known as Voronoi tesselation [28], constituted by convex polygons, each one composed of the points closer to a Delaunay vertex (its center) than to any other. Voronoi polygons naturally appear in many situations in nature, such as the hottest parts of the sun, ice melting, and even on giraffe skin. Apart from these relations connected with natural science, not directly related to the numerical solution of PDEs, more relevant for us is the use of Voronoi grids for the discretization by finite volumes for oil reservoir simulations, oceanographic models, and other complex problems in science and engineering which are derived from conservation laws. Finite volume methods are discretizations locally preserving some conservation properties. A first approach to this kind of methods was introduced for one-dimensional problems by Samarskii in 1960,see [24], as a finite difference method called balance method (or integro-interpolation method). For many years finite volume methods on rectangular and triangular grids were used in heat transfer and computational fluid mechanics [6, 20, 22]. From the theoretical point of view, we refer to the reader to Eymard et al. [9]. One of the most important aspects in the numerical solution of partial differential equations is the efficient solution of the corresponding large systems of equations arising from their discretization. Multigrid methods [3, 11, 27] are among the most powerful techniques for solving such type of algebraic systems, and they have become very popular among the scientific community. Basically, multigrid methods are based on the property of a strong smoothing effect on the error of many iterative methods, together with the fact that an smooth function can be well represented on coarser grids, where its approximation is less expensive. The multigrid performance strongly depends on the choice of the components of the algorithm, which are problem dependent, and it is crucial an harmonic interplay between the smoothing and the coarse-grid correction. There are two basic approaches to multigrid solvers: geometric and algebraic 2 multigrid. Whereas in geometric multigrid a hierarchy of grids must be proposed, in algebraic multigrid no information is used concerning the grid on which the governing PDE is discretized. The focus on this work is to consider geometric multigrid methods. Not many authors have applied geometric multigrid methods on cell-centered discretizations. Most of these works have been done in rectangular grids. Historically, the pioneer work is due to Wesseling, [31], in which a multigrid method for interface problems was constructed to simulate oil reservoir problems. This work started a chain reaction of papers focused on this subject, see [12, 13, 14, 15, 29, 30, 32]. The W-cycle convergence of these multigrid methods, in the case of natural injection as prolongation, was theoretically analyzed by Bramble et al. [2], and in the case of V-cycle, the convergence was proved by Kwak et al. in [17] using certain weighted prolongation operators. By other hand, on triangular grids, Kwak et al., [18, 19], proposed a new multigrid method, extending their previous works. For an irregular domain, it is very common to apply regular refinement to an unstructured input grid. In this way, a hierarchy of locally structured grids is generated. To perform this refinement, each triangle is divided into four congruent ones by connecting the midpoints of their edges. In this way, a hierarchy of grids is obtained, where transfer operators between two consecutive levels can be defined. These grids provide a suitable framework for the implementation of a geometric multigrid algorithm, permitting the use of stencil-based data structures, see [1], being necessary only a few stencils to represent the discrete operator, reducing drastically the memory required. In this work, very simple local inter-grid transfer operators have been chosen to make easier the communication between different input blocks. For this reason, powerful smoothers have to be designed. Different smoothers are used on each input block of the initial unstructured grid, depending on its geometry. Notice that these smoothers couldn’t be implemented on a pure unstructured grid. We are speaking about, for example, the line-type smoothers, very necessary in the case of anisotropic problems. Regarding particular applications, the diffusion equation with discontinuous coefficients is widely used in computational fluid dynamics. However, it is well-known that when standard multigrid is applied for solving equations with highly varying discontinuous coefficients, a deterioration in the convergence of the method can be obtained, and even divergence can be observed. In the case of cell-centered approximations, the main idea to circumvent the problem is to use Galerkin approximation [32] on coarse grids. The disadvantage usually observed when using this approach, is that the stencils of the coarsegrid operators are often larger than the corresponding fine-grid stencil, what is problematic especially in three dimensions. However, note that the use of simple transfer operators preserves four-point stencils in the case of cell-centered schemes on triangular grids. In the case of composite materials their components have nearly constant diffusivity, but vary by several orders of magnitude. In these cases, it is quite common to idealize the diffusivity by a piecewise constant function. This makes suitable the use of semi-structured grids to deal with these problems. 3 The rest of the work is organized as follows: Chapter 2 presents the basic principles of Multigrid. Chapter 3 is devoted to present the cell-centered discretization, deriving the corresponding equations which define the discrete operator for both unstructured and structured triangular grids. A multigrid method for this kind of scheme is proposed in Chapter 4. In particular, special smoothers are designed to overcome troubles arising from some difficult grid-geometries. In Chapter 5, the good performance of the proposed multigrid methods is demonstrated on structured triangular grids, comparing the behavior of the considered smoothers for different grid-geometries. In Chapter 6, some test examples are presented to show the applicability of the strategies followed throughout all the work. Finally, in Chapter 7, we focus the previous work in the resolution of diffusion equation with piece-wise discontinuous coefficients. 4 Chapter 2 Introduction to Multigrid Multigrid methods, are a family of iterative methods based on the property of some classical iterative algorithms, which can annihilate the high frequency components of the error after few iterations, which leads to an smoothed error, see Figure 2.1. Figure 2.1: Representation of the error after: 0,5,10 and 50 iterations of a classical iterative algorithm. Considering a linear problem, Au =b, (2.1) where Ais an n×nmatrix, uthe unknowns vector and bthe right hand side. Let us consider the following splitting notation: A=M−N, (2.2) where Mis an easily invertible matrix. By substituting in 2.1 we obtain (M− N)u=bor equivalently, Mu =Nu+b. Now, by denoting umas the approximate solution after the m-iteration, the general scheme for an iterative algorithm is: um+1 =M−1Num+M−1b(2.3) 5 If we call ueto the exact solution, we can define the error in the iteration m, as em=ue−um,and the residual rm=b−Aum.We can combine these expressions with (2.1) in order to get a system of equations for the error: rm=b−Aue+Aem=Aem,(2.4) therefore, we can solve a system for the error in order to get a better approximation of the solution. If we now combine the smoothing property, that allows to represent the error in a coarser grid without data loss, with the possibility to solve an equation for the error, we can define the following strategy: 1◦Smooth. 2◦Represent the error in a coarser grid. 3◦Solve the equation exactly. 4◦Use the error calculated to correct the approximation, um. 5◦Iterate until obtaining the desired precision. The above explained method is called two-grids method. Its efficiency is not a high improvement, however, its convergence rate is constant not matter the number of nodes, this is also known as h-independence rate. Nevertheless, two-grids method can be improved if, in the third step, we approximate the analytical solution by using another two-grids method and so on, leading to the so called multigrid method, which keeps the h-independence rate, and its numerical cost is about O(n×n). Regarding the second and the forth step we can split the inter-grid transfer operators in: •Restriction: Introduce, the data from a fine grid to a coarse one, it can be done by simple injection, just introduce in the coarse grid the value of the nodes that relies in the same coordinates, or by an average of nodes. •Prolongation: Interpolate the solution of the coarse grid to the nodes of a finer one. The simplest, the adjoint of the simple injection, is to introduce the data from the coarse grid directly to a cluster of nodes. One can realize that the performance of the multigrid method relies on the smoother, coarse grid transfer method and in the interpolation from a coarser grid to a finer one. The intercourse between the smoother and the inter-grid transfer operators are the key in multigrid. In fact there are two families of multigrid methods, depending on their intercourse: •Geometric multigrid: The coarsening is fixed, usually the coarser grid is twice the finer, and the smoother must annihilate the half of the frequency errors. These methods are fast, but the geometry must be simple. •Algebraic multigrid: The smoother is fixed, and the coarsening must adapt itself in order to compensate what the smoother cannot smooth. These family fits very easily complex geometries, but they are slower than the geometric multigrid. 6 Figure 2.2: Relation, along the grids, between nodes for a simple injection and its adjoint. Regarding the smoothers, the most common smoothers are the Jacobi iterative method with a relaxation parameter, Gauss-Seidel algorithm and a modification of the Gauss-Seidel, in which you apply a Gauss-Seidel step but in a chessboard manner, i.e: first the odd nodes and next the even ones, see Figure 2.3. Figure 2.3: Red-Black smoother on a triangular grid. There are three kind of multigrid methods depending on the effort to approximate the error on coarser grids, one get three different multigrid methods, V-cycle, W-cycle and F-cycle, see Figure 2.4. V-cycle is the fastest, but the weakest, F-cycle is usually preferred over W-cycle since it achieve similar re- 7 Figure 3.6: Notation for neighboring Voronoi points on an structured grid, characterized by angles αand β. and the lengths of the sides of Tare l2=l, l1=lsin β sin(α+β),and l3=lsin α sin(α+β), after simple calculations, we finally obtain the stencils: Ldu h=2 sin(α+β) l2sin αsin β  0 tan(α+β) 0 −tan α−tan β0 0 0 0  ,(3.16) Ldd h=2 sin(α+β) l2sin αsin β  0 0 0 0 tan(α+β)−tan α−tan β0 0 0 0  .(3.17) As previously commented, for an arbitrary up-oriented Voronoi grid-point xu k1,k2 in Vu hthe stencil would be the “mirror image stencil” of (3.16)-(3.17), that is, Luu h=2 sin(α+β) l2sin αsin β  0 0 0 0 tan(α+β)−tan α−tan β0 0 0 0  ,(3.18) Lud h=2 sin(α+β) l2sin αsin β  0 0 0 0−tan β−tan α 0 tan(α+β) 0  .(3.19) Notice that depending on the angles characterizing the grid, some strong connections appear between unknowns and this fact makes that we will have to take care in the design of the smoothers in a geometric multigrid method, as will be discussed in next chapter. 14 Chapter 4 Multigrid method The performance of geometric multigrid methods is strongly dependent on the choice of adequate components to the considered problem. The main components are the smoother Sh,inter-grid transfer operators: restriction I2h hand prolongation Ih 2h,and the coarse-grid operator L2h.These components have to be chosen so that they efficiently interplay with each other in order to obtain a good connection between the relaxation and the coarse-grid correction. In this chapter the proposed cell-centered multigrid algorithm is described. All the attention is focused in the detailed explanation of the considered smoothers and the special features appearing due to the cell-centered character of the discretization. Although the presentation of such components is done on a regular structured grid, our purpose is to apply the proposed multigrid method in the framework of semi-structured grids. Therefore, the choice of the corresponding components is done also with a view to this application. In this case, we will use a block-wise multigrid algorithm, where each triangle of the coarsest grid is treated as a different block with regard to the smoothing process. This blockwise strategy is suitable thanks to the possibility of choosing different smoothers for triangles having different geometries, thus resulting in an improvement of the characteristics of our algorithm. Besides, we will have to take care in the communication among the triangles of the coarsest triangulation. Next, we are going to describe the components of the algorithm that we are going to consider throughout all this work. 4.1 Coarse-grid correction In the application of geometric multigrid, a hierarchy of grids is needed in order to accelerate the convergence of the smoother, by using solutions obtained on the coarser meshes as corrections. As previously commented, in order to obtain such hierarchy of grids, we divide the initial triangles into four congruent ones by connecting the midpoints of the edges, and so forth until the mesh has the desired fine scale to approximate the solution of the problem. When vertex-centered discretizations are considered on triangular grids, grid points laying on coarser grids also belong to the finer grids, giving rise to a 15 so-called nested hierarchy of grids. However, when the considered cell-centered discretizations are used, it is worth to note that except in the case of equilateral triangles, the grid hierarchy results to be non-nested, see Figure 4.1. Figure 4.1: Nested (left) and non-nested (right) hierarchies. This makes the interplay between smoothing and coarse-grid correction specially difficult, requiring the design of new smoothers or new inter-grid transfer operators. With a view to the application of the proposed multigrid method on semi-structured grids, our proposal is to consider very simple inter-grid operators, since this will facilitate the communication between the blocks composing these grids. In particular, injection operator is considered as the prolongation, and its adjoint is chosen as the restriction, resulting in the fact that only the four fine-grid points surrounding a coarse-grid point contribute to its restriction. More concretely, the considered restriction operator, I2h h,is given in the following way I2h h=  (I2h h)uu (I2h h)ud (I2h h)du (I2h h)dd  ,with                  (I2h h)uu =  0 0 0 0 1/4 0 0 0 0  ,(I2h h)ud =  0 0 0 0 1/4 1/4 0 1/4 0  , (I2h h)du =  0 1/4 0 1/4 1/4 0 0 0 0  ,(I2h h)dd =  0 0 0 0 1/4 0 0 0 0  , as shown in Figure 7.9, and the corresponding prolongation fulfils Ih 2h= 4 I2h h. The choice of these inter-grid transfer operators leads us to make an effort in the Figure 4.2: Restriction operator. smoothing process. Then, we must design efficient smoothers capable of taking charge of the remaining components of the error, which cannot be eliminated by the coarse-grid correction part of the algorithm. Actually, the design of suitable smoothers for cell-centered grids is a challenge in this context. 16 4.2 Smoothers The smoother usually plays an important role in multigrid algorithms, mainly in the geometric approach, resulting the choice of a suitable smoother an important feature for the efficiency of these methods. Moreover, as previously commented, in the framework we are working with, this choice is even more relevant. Due to the general observation that errors become smooth if strongly connected unknowns are collectively updated, appropriate smoothers have been designed depending on the magnitude of the coefficients of the stencils, given in (3.16)- (3.19). The following smoothers have been considered and tested in order to fulfil the previous requirement. 4.2.1 Jacobi smoother. For almost equilateral triangles, the magnitude of all the entries of the stencils is similar, and therefore a point-wise smoother is enough to satisfactorily reduce the components of the error. The easiest smoother to perform is a Jacobi type smoother, which consists of computing the approximation of each unknown, by using non-updated values of the rest of the unknowns. Notice that this implies that it makes not difference the order in which the grid points are visited, what makes Jacobi scheme well-suited for parallel processing. However, for difficult problems, usually this smoother does not give enough satisfactory results, and some variants have to be considered. As is the case for Jacobi smoother presented here, some standard smoothers are based on a decomposition on the positive and negative parts of the operator, which correspond to the updated an non-updated unknowns before the current step. Here we are going to present the corresponding decomposition for Jacobi smoother. In order to do this, only the positive part of the operator is displayed. Taking into account that only the diagonal blocks of the operator contribute in this positive part, it holds (Luu h)+=Luu h,(Ldd h)+=Ldd h.(4.1) 4.2.2 Red-Black smoother. Due to the fact that unknowns related to up or down-oriented triangles have no direct connection with each other, it seems natural to simultaneously update all unknowns associated with equally oriented triangles, giving rise to a pattern relaxation scheme. Since two different types of grid-points are distinguished, a two-color relaxation process, called here Red-Black smoother, is considered. More concretely, one iteration of this relaxation scheme consists of two partial steps. In the first one, unknowns corresponding to up-oriented triangles are updated, and in the second step those associated with the down-oriented triangles are relaxed by using the updated values. Thus, the complete smoothing operator Shis given by the composition of two partial step operators, Su hand Sd h,which correspond to apply a Jacobi step on each type of grid-points, that is, Sh=Sd hSu h.These partial step operators are characterized by a different decomposition than the previous Jacobi over all the grid-points, in the way that for Su h,for example, the positive parts of the scalar operators are: (Luu h)+=Luu h,(Ldd h)+=Ih,(4.2) 17 and if Sd his considered, the identity operator will correspond to (Luu h)+,and (Ldd h)+=Ldd h. 4.2.3 Diamond smoothers. For almost rectangular triangles, a strong connection between only two nodes involved in the stencil appears, due to the anisotropy of the Voronoi mesh. Therefore, since the common lore claims that smoothing must be done in the direction of the strong connection, in this case, both unknowns will have to be simultaneously relaxed. These unknowns are associated with the closest Voronoi points corresponding to different-oriented triangles, as seen in Figure 4.3a). Therefore, a small (2 ×2)-system must be solved for each of these pairs of unknowns. Different orderings can be chosen to visit these blocks. We have chosen the lexicographic one, but of course, many orderings are possible. In triangular grids, three different diamond smoothers, associated with the three edges of a triangle, can be defined. If a triangle characterized by angles αand β is considered, we can assign a different color to each of its vertex, in the way that, for example, black color is associated with the vertex corresponding to angle α, green color with that vertex associated with β, and red color corresponds to vertex of angle 180o−(α+β).In this manner, each diamond smoother can be named with the color corresponding to the vertex opposite to its associated edge. Following this rule, diamond smoother appearing in Figure 4.3a) is called green-diamond smoother. a) b) Figure 4.3: a) Green-diamond smoother and b) Red-wormy smoother. This smoother is based on a decomposition of positive and negative parts of the operator. Although three different diamond smoothers can be considered, the corresponding decomposition can be done analogously. Then, in order to do this description, we consider the green-diamond smoother, which can be seen in Figure 4.3a). The positive parts of the involved scalar operators are given as 18 follows: (Luu h)+=Luu h,(Lud h)+=  000 0sud 0,00 0sud 0,−10  , (Ldu h)+=  0 0 0 sdu −1,0sdu 0,00 0 0 0  ,(Ldd h)+=Ldd h. (4.3) 4.2.4 Wormy smoothers. When a very small angle characterizes the triangulation, the strong coupling appears between the Voronoi points associated with the up and down-oriented triangles connected in the direction of the edge opposite to the vertex corresponding to this small angle. Therefore, all those points are simultaneously updated, see Figure 4.3b), and a tridiagonal system must be solved for each of these “wormy-lines”. For this reason, this smoother will be called wormysmoother. Similarly to the previous case of the diamond smoother, in triangular grids, three different wormy-smoothers can be defined associated with the three edges of a triangle. Again, each wormy-smoother can be named with the color corresponding to the vertex opposite to its associated edge. Following this criterion, wormy-smoother appearing in Figure 4.3b) is called red-wormy smoother. Notice that these smoothers can be performed by visiting the “lines” from vertex to edge or, on the contrary, from edge to vertex, and this latest is the chosen option. From the description of this smoother, we can obtain the decomposition of the discrete operator which gives rise to wormy-smoother. Analogously to the diamond smoother, three different wormy smoothers can be defined, whose decomposition can be obtained in a similar way. In order to present the corresponding decomposition of the operator, the red-wormy smoother, displayed in Figure 4.3b), is considered. In this way, the positive parts of the scalar operators result as follows: (Luu h)+=Luu h,(Lud h)+=  0 0 0 0sud 0,0sud 1,0 0sud 0,−10  , (Ldu h)+=  0 0 0 sdu −1,0sdu 0,00 0 0 0  ,(Ldd h)+=Ldd h. (4.4) Concluding, we can say that each of these wormy-smoothers will be suitable when the angle corresponding to the vertex of its color is small, and in this way, any possible triangulation will have associated a wormy-smoother giving a satisfactory convergence factor. 19 Chapter 5 Results of the proposed multigrid method on structured grids In order to see the suitability of the introduced smoothers in the design of an efficient multigrid algorithm, we present some experiments comparing their behavior in regular structured grids. In particular, model problem (3.1)-(3.2) is solved in the structured grid arising from the regular refinement of a triangular domain, characterized by two of their angles. In all the presented experiments, an F(2,2)−cycle is considered. F−cycle is preferred to V−cycle due to the poor chosen inter-grid transfer operators. We begin by considering an equilateral triangular domain. As any anisotropy arises from the grid resulting of a regular refinement, it could be seen natural to think in applying a simple point-wise smoother, like Jacobi or lexicographic Gauss-Seidel, to this kind of triangulations. However, as previously commented, for this situation a pattern relaxation scheme could be more appropriately. Thus, we are going to present some convergence results on a regular equilateral grid, comparing the behavior of multigrid by considering: undamped Jacobi smoother, Red-Black smoother, ω−Red-Black smoother (with ω= 1.15) and diamond smoother. In Figure 5.1a), we show the history of the convergence on a grid obtained after eight refinement levels, and by considering as stopping criterion to reduce the maximum residual until 10−8.First of all, a rather surprising observation could be concluded from this figure: the performance of undamped Jacobi appears to be a satisfactory choice as smoother for cell-centered discretizations on triangular grids (as also seen for other type of discretizations on triangular grids, [10], and in the context of full-multigrid on rectangular grids, [23]), despite the well-known lack of smoothing property of this iterative scheme. Notwithstanding this unusual behavior, the obtained Jacobi results are largely improved by Red-Black smoother and diamond smoother. At the same time, the convergence factors provided by both smoothers are enhanced by the Red-Black smoother with relaxation parameter ω= 1.15,which has been obtained by experimental tests. This improving effect was pointed out in [21] for cell-centered discretizations on rectangular grids, where it was validated by a local Fourier analysis. Moreover, the good behavior of the multigrid based 20 a) b) Figure 5.1: a) Comparison of smoothers on an equilateral triangular grid. b) History of the convergence for different numbers of refinement levels by using ω−Red-Black smoother. on the ω−Red-Black smoother is confirmed in Figure 5.1b), where its robustness with regard to the discretization parameter is shown. In this figure, the history of the convergence of the method is displayed for different numbers of refinement levels, resulting to be independent. Therefore, we conclude from this experiment that ω−Red-Black smoother seems to be a good choice for almost equilateral triangulations. However, this good behavior deteriorates very quickly when the shape of the triangle tends to be rectangular or is characterized by a very small angle. In the case of almost rectangular triangles, point-wise smoothers are not suitable anymore due to the anisotropy of the Voronoi mesh. By other hand, diamond smoother results in a very efficient smoother when this kind of grids are considered. As an example, a triangular domain characterized by angles α= 45oand β= 85ois fixed. In Figure 5.2a), the history of the multigrid convergence by using different smoothers is displayed. More concretely, Red- Black, ω−Red-Black and diamond smoothers are used in this comparison. In all cases, the finest grid results by applying eight refinement levels to the initial triangular domain. a) b) Figure 5.2: a) Comparison of smoothers on an almost rectangular triangular grid. b) History of the convergence for different numbers of refinement levels by using diamond smoother. 21 As we can observe, very poor rates are obtained when both Red-Black smoothers are considered, whereas the convergence factors provided by the new diamond smoother are very satisfactory, achieving the convergence in only eleven cycles. Besides, in Figure 5.2b), where the history of the convergence is shown for different numbers of refinement levels, the robustness of this smoother with respect to the space discretization parameter is demonstrated. Although convergence factors provided by diamond smoother are very satisfactory for many grid configurations, when a triangulation characterized by a very small angle is used, this smoother gives rise to poor rates. This behavior can be seen in Figure 5.3, where asymptotic convergence factors of the diamond smoother based multigrid are shown for a wide range of pairs of angles characterizing the grid. Figure 5.3: Experimentally computed convergence factors for the diamond smoother based multigrid and four smoothing steps, for different triangles in function of two of their angles. To overcome these troubles appearing when the primal mesh is anisotropic, wormy-smoother in the direction of the anisotropy is a suitable smoother, largely improving the convergence factors provided by the rest of point-wise or blockwise smoothers. To validate this statement, we are going to compare the multigrid convergence by using each one of the smoothers proposed in this work, when an isosceles triangle with a small angle of 10ois considered as domain of our problem. With this purpose, in Figure 5.4a), the multigrid convergence provided by using ω−Red-Black, diamond and wormy smoothers is depicted for eight refinement levels. From this picture, it is clear that wormy-smoother is the best choice for this type of triangulations. Moreover, an h−independent convergence is also shown in Figure 5.4b). Concluding, from the results presented in this chapter, it seems that a reasonable strategy to follow would be to apply the point-wise ω−Red-Black smoother for almost equilateral triangles, the collective diamond smoother for almost rectangular triangles, and finally the appropriate block collective wormy smoother 22 a) b) Figure 5.4: a) Comparison of smoothers on an isosceles triangular grid with smallest angle 10o.b) History of the convergence for different numbers of refinement levels by using wormy-smoother. when triangulations with a small angle appear. 23 We propose an easy algorithm to calculate Lk−1using only stencils. This method must be applied once per stencil operator. Firstly, the notation for this algorithm needs to be modified, nevertheless using the following conversion of the previously explained stencils, its application is straightforward: Stencil =  000 Lud h(0,0) Luu h(0,0) Lud h(1,0) 0Lud h(0,1) 0  .(7.8) Secondly, we have to define an stencil of stencils: St(j, i) =   0 0 0 SE(jj, ii)SC(jj, ii)SW(jj, ii) 0SN(jj, ii)|SS(jj, ii) 0  ,(7.9) Where the subfix is their position in cardinal coordinates, being SCthe center, see Figure 7.2. Figure 7.2: Stencils notation to be used in the Stencil of Stencils All the stencils coordinates must be defined from −1 to 1,being the center (0,0), and the κmust be included. The stencil of stencil, Stmust be a 4 dimensional array in which the two first coordinates are the position of the stencil, and the last two the coordinates the position entries of that stencil. The algorithm is as follows: Program Galerkin_Operator StencilC = 0 n=0 center = 1 !Compute A_coarse doii=0to1 do jj = 1 to -1 add -1 doi=0to1 do j = 1 to -1 add -1 n=n+1 if((ii=0andjj=0)or(i=0andj=0)or (n = center)) then StencilC(0,0) = StencilC(0,0) + + S_t(j,i,-jj,ii) * R(j,i) 30 else StencilC(jj,ii) = StencilC(jj,ii) + + S_t(-j,i,jj,ii) * R(-j,i) end do end do n=0 center = center + 1 end do end do !Store stencil StencilC = StencilC * Cf Where Cf, is an scale constant to compensate the fact that the local prolongation and restriction do not satisfy the following inequation [32]: mp+mr> m. (7.10) In this equation, mpand mrare the order of the prolongation and the restriction, respectively, that is, the highest degree plus one of polynomials that are interpolated exactly, and mis the order of the differential equation. When the simple injection is selected, the Galerkin operator is not consistent and it must substituted by the following expression [33] Lk−1=1 2Ik−1 kLkIk k−1, k = 1, . . . , f (7.11) this modified version provides better multigrid convergence rate results over the original Galerkin operator. 7.2 Numerical experiments In this section we are going to present two numerical experiments to demonstrate the efficiency of the proposed block-wise multigrid algorithm based on Galerkin approach. In the first experiment, problem (7.1) is solved on the unit square with two different distributions of diffusion coefficients. This is a typical benchmark problem on this study area. In the second one, the same problem is solved on a more complex domain, simulating a more real situation. In both experiments the proposed multigrid algorithm is applied by using an F-cycle with two pre- and two post-smoothing steps. 7.2.1 Diffusion problem on the unit square with various diffusion coefficients In the first experiment, we are going to solve problem (7.1) on the unit square with different diffusion coefficients defined piecewise on two subdomains of different shapes, as shown in Figures 7.3 and 7.4. This consists of a benchmark problem taken from the literature [32]. More concretely, in the first test case, the inner subdomain has a rhombus shape with a 0.5 side length, see Figure 7.3, whereas the second one consists of an hexadecagon with diameter 0.5 simulating a circle, see Figure 7.4. The right-hand side is defined as f=xy, the Dirichlet boundary conditions are given by g=x2+y2,and the diffusion coefficient values 31 are 0.333 ×105for the internal subdomains and 2 for the rest of the square, see Figures 7.3 and 7.4. In the same figures the corresponding coarsest grids are also represented. a) b) Figure 7.3: a) Coarsest unstructured mesh for the first test case, and distribution of diffusion coefficients: 0.333×105at the white region and 2 at the shaded part. b) Different smoothers for the triangles of the coarsest grid: white is Red-Black smoother, diamond smoother is represented by light-grey, and wormy smoother by dark-grey. a) b) Figure 7.4: a) Coarsest unstructured mesh for the second test case, and distribution of diffusion coefficients: 0.333 ×105at the white region and 2 at the shaded part. b) Different smoothers for the triangles of the coarsest grid: white is Red-Black smoother, diamond smoother is represented by light-grey, and wormy smoother by dark-grey. For these examples, Red-Black, Wormy and Diamond smoothers have been considered, along with some extra-relaxation process on Voronoi nodes close to the internal boundaries of the initial coarsest grid. The proposed block-wise multigrid method has been applied to solve both test cases. Red-black, wormy and diamond smoothers have been used for different triangles of the coarsest grid, as shown in Figures 7.3 7.4. Regarding the obtained multigrid convergence, in Table 7.2.1, the number of iterations necessary to reduce the initial residual in a factor of 10−10 are shown for both test cases, together with the asymptotic convergence rates and the number of unknowns for each number of refinement levels. We observe an h-independent convergence for both problems, and although these results are slightly worse than those obtained in the case of constant diffusion coefficients, as expected, 32 Diamond Circle Levels Unknowns Cycles ρhCycles ρh 4 6912 8 0.15 8 0.13 5 27648 8 0.18 9 0.17 6 110592 8 0.22 9 0.19 7 442368 8 0.24 9 0.20 8 1769472 9 0.25 9 0.21 Table 7.1: Number of iterations to reduce the residual in a factor of 10−10,and asymptotic convergence rates for both test cases. the method shows a very satisfactory convergence. On the other hand, when direct discretization is used on coarse grids, a very poor convergence rate is obtained. 7.2.2 Diffusion problem on a composite material domain In the second experiment, problem (7.1) is solved on a rectangular domain composed of two different materials with different diffusion coefficients: 1 and 0.001, as we can see in Figure 7.5a). The considered coarsest grid is also shown in the same figure, and also we can observe that it is composed of triangles with very disparate geometries. For this reason different smoothers are considered for the triangles of the coarsest triangulation. In particular the smoothers chosen for these triangles are shown in Figure 7.5b). In this way, the block-wise multigrid previously proposed is used for solving the problem. We want to comment that an extra-relaxation on Voronoi nodes close to the internal boundaries of the initial coarsest grid has been necessary. a) b) Figure 7.5: a) Diffusion coefficients for the second experiment, white color represents κ= 0,001 and grey κ= 1. b) Different used smoothers: white is red-black, diamond smoother is represented by light-grey and wormy smoother by dark-grey. Now we are going to compare the behavior of the multigrid algorithm by considering both, direct discretization on coarse grids and the Galerkin approach. For this purpose, in Figure 7.6, the history of the convergence of the method, for two different refinement levels, to reach a final maximum residual below 10−7, is displayed. In this case, the method based on direct discretization leads to divergence, while that based on Galerkin approach yields very satisfactory and robust results. Moreover, we observe that the convergence is independent of the discretization parameter, and with only twelve/thirteen cycles the residual 33 reaches the desired value. Note, that in this experiment the use of Galerkin coarse-grid operator is mandatory. Figure 7.6: Comparison between direct discretization and Galerkin for two different refinement levels. 34 Chapter 8 Conclusions In this work, the design of efficient multigrid methods for cell-centered finite volume schemes on semi-structured triangular grids has been the main focus. Due to the cell-centered character of the discretization and because of the application of the proposed strategy on semi-structured grids, very simple inter-grid transfer operators have been used in the design of these methods, leading to the requirement of stronger smoothers. Thus, appropriate novel smoothers are proposed depending on the geometry of the grid. Moreover, due to the semistructured nature of the grid, different smoothers can be used on each structured region, giving rise to a block-wise multigrid method. The global behavior relies on its components on each block. Furthermore, the good behavior of the proposed multigrid method has been illustrated by some numerical experiments. A fast and h-independent convergence has been obtained, concluding that the adopted strategy yields very efficient solvers on relatively complex domains for cell-centered discretizations. Also, the proposed Galerkin operator gives good convergence rates even for high diffusion coefficient jumps, concluding that this methodology is efficient also for heterogeneous materials, which, in first instance are prohibitive. Moreover, the algorithm proposed is fast an efficient as it only use the collide stencils to get the Galerkin operator, making the user of this operator more plausible than having to multiply three matrices. 35 Bibliography [1] B. Bergen, T. Gradl, F. H¨ ulsemann, U. Ruede,A massively parallel multigrid method for finite elements, Comput. Sci. Eng., 8: 56–62, 2006. [2] J. Bramble, R. Ewing, J. Pasciak, J. Shen,The analysis of multigrid algorithms for cell centered finite difference methods, Adv. Comput. Math., 5: 15-29, 1996. [3] A. Brandt,Multi–level adaptive solutions to boundary-value problems, Math. Comput., 31: 333–390, 1977. [4] J. Brandts, S. Korotov, M. Kˇ r´ ıˇ zek,The discrete maximum principle for linear simplicial finite element approximations of a reaction–diffusion problem, Linear Algebra Appl., 429: 2344–2357, 2008. [5] P.G. Ciarlet,The finite element method for elliptic problems, Series “Studies in Mathematics and its Applications”, North-Holland, Amsterdam, 1978. [6] P. Colella, E.G. Puckett,Modern Numerical Methods for Fluid Flow, Lecture Notes, Department of Mechanical Engineering, University of California, Berkeley, CA, 1994. [7] H. Erten, A. ¨ Ung¨ or,Quality Triangulations with Locally Optimal Steiner Points, SIAM J. Sci. Comput., 31: 2103–2130, 2009. [8] H. Erten, A. ¨ Ung¨ or,Computing Triangulations without Small and Large Angles, International Symposium on Voronoi Diagrams (ISVD) 2009. [9] R. Eymard, T. Gallou¨ et, R. Herbin,Finite volume methods, in: Handbook of Numerical Analysis, Vol.7, North Holland, Amsterdam, 2000, pp. 713–1020. [10] F.J. Gaspar, J.L. Gracia, F.J. Lisbona,Fourier analysis for multigrid methods on triangular grids, SIAM J. Sci. Comput., 31: 2081–2102, 2009. [11] W. Hackbusch,Multi–grid methods and applications, Springer, Berlin, 1985. [12] M. Khalil,Analysis of linear multigrid methods for elliptic differential equations with discontinuous and anisotropic coefficients. PhD thesis, Technical University Delft, The Netherlands, 1989. 36 [13] M. Khalil, P. Wesseling,A cell-centered multigrid method for threedimensional anisotropic-diffusion and interface problems, in: Mandel, J., McCormick, S.F. (eds.), Preliminary Proc. of the 4th Copper Mountain Conference on Multigrid Methods, Vol. 3, pp. 99-117, Denver, 1989. Computational Mathematics Group, Univ. of Colorado. [14] M. Khalil, P. Wesseling,Vertex-centered and cell-centered multigrid for interface problems, in: Mandel, J., McCormick, S.F. (eds.), Preliminary Proc. of the 4th Copper Mountain Conference on Multigrid Methods, Vol. 3, pp. 61-97, Denver, 1989. Computational Mathematics Group, Univ. of Colorado. [15] M. Khalil, P. Wesseling,Vertex–centered and cell–centered multigrid for interface problems, J. Comput. Phys., 98: 1-20, 1992. [16] J. Kar´ atson, S. Korotov,Discrete maximum principles for finite element solutions of nonlinear elliptic problems with mixed boundary conditions, Numer. Math. 99: 669-698, 2005. [17] D.Y. Kwak,V–cycle multigrid for cell centered finite differences, SIAM J. Sci. Comput., 21: 552-564, 1999. [18] D.Y. Kwak, J.S. Lee,Multigrid algorithm for the cell centered finite difference method II: Discontinuous coefficient case, Numer. Meth. Partial Differential Eqs., 20(5): 723-741, 2004. [19] D.Y. Kwak, J.S. Lee,Comparison of V–cycle Multigrid Method for Cell– Centered Finite Difference on Triangular Meshes, Numer. Meth. Partial Differential Eqs., 22: 1080–1089, 2006. [20] R.J. Leveque,Finite volume methods for hyperbolic problems, Cambridge Texts in Applied Mathematics. Cambridge: Cambridge University Press. xix, 558, 2002. [21] M. Mohr, R. Wienands,Cell–centered multigrid revisited, Comput. Visual. Sci., 7: 129–140, 2004. [22] S.V. Patankar, D.B. Spalding,Heat and mass transfer in boundary layers, Morgan–Grampian, London, 1967. [23] C. Rodrigo, F.J. Gaspar, C.W. Oosterlee, I. Yavneh,Accuracy Measures and Fourier Analysis for the Full Multigrid Algorithm, SIAM J. Sci. Comput., 32: 3108–3129, 2010. [24] A.A. Samarskii,Parabolic equations and difference methods for their solution, in: Proc. of All–Union Conference on Differential Equations, 1958, Publishers of Armenian Ac.Sci., 1960, pp 148–160. [25] J.R. Shewchuk,Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator, in “Applied Computational Geometry: Towards Geometric Engineering” (Ming C. Lin and Dinesh Manocha, editors), volume 1148 of Lecture Notes in Computer Science, pages 203–222, Springer- Verlag, Berlin, May 1996. (From the First ACM Workshop on Applied Computational Geometry.) 37 [26] J.R. Shewchuk,Delaunay Refinement Algorithms for Triangular Mesh Generation, Comput. Geom., 22: 21–74, 2002. [27] U. Trottenberg, C.W. Oosterlee, A. Sch¨ uller,Multigrid, Academic Press, New York, 2001. [28] C. Voronoi,Nouvelles applications des param`etres continus `a la th´eorie des forms quadratures, J. Reine Angew Math., 134: 198–287, 1908. [29] P. Wesseling,Cell–centered multigrid for interface problems, J. Comput. Phys. 79: 85-91, 1988. [30] P. Wesseling,Cell–centered multigrid for interface problems, in: Mc– Cormick, S.F. (ed.), Multigrid Methods: Theory, Applications, and Supercomputing, Vol. 110 of Lecture Notes in Pure and Applied Mathematics. New York: Marcel Dekker 1988, pp. 631-641. [31] P. Wesseling,A multigrid method for elliptic equations with a discontinuous coefficient, in: van der Burgh, A.H.P., Mattheij, R.M.M. (eds.), Proceedings of the First International Conference on Industrial and Applied Mathematics (ICIAM 87). Philadelphia: SIAM 1988, pp. 173–183. [32] P. Wesseling,An introduction to multigrid methods, John Wiley, Chichester, UK, 1992. [33] D.J. Mavriplis, V.Ventakatrishnan,Agglomeration multigrid for viscous turbulent flows, ICASE, 1994 38