Repositorio Institucional de Documentos

2013 128 Pablo Salinas Cortés Semi-Structured multigrid methods on Voronoi meshes to the resolution of the Darcy-OberbeckBoussinesq model Departamento Director/es Ingeniería Mecánica Gaspar Lorenz, Francisco José Lisbona Corté, Francisco Javier Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es Pablo Salinas Cortés SEMI-STRUCTURED MULTIGRID METHODS ON VORONOI MESHES TO THE RESOLUTION OF THE DARCY-OBERBECK-BOUSSINESQ MODEL Director/es Ingeniería Mecánica Gaspar Lorenz, Francisco José Lisbona Corté, Francisco Javier Tesis Doctoral Autor 2013 Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Semi-Structured multigrid methods on Voronoi meshes to the resolution of the Darcy-Oberbeck-Boussinesq model Pablo Salinas Cort´es Applied mathematics University of Zaragoza A thesis submitted for the degree of Doctor (PhD) of Computational Mechanics 2013 June ii To my beloved Cristina Gargallo Acknowledgements Regarding the mathematical team, firstly I would like to thank Francisco Gaspar for offering me the possibility of doing a PhD, and for his help during these years whenever I needed it. Also, I want to expose my gratitude to Francisco Lisbona since he always explained me any problem I faced, he is like the encyclopaedia of the group. And last but not least, I want to give many thanks to Carmen Rodrigo because she always was there for quick doubts, except when she was in the other side of the world. I want to emphasize my gratitude for their help correcting this thesis. Special thanks to Kees Oosterlee for hosting me in Amsterdam and showing me Delft during my staying in the Netherlands, where I had an incredible experience. I would like to express my gratefulness to my family, Javier, Margarita and David, for their support, patience and for cheering me on during this project despite at the beginning it seemed very difficult to fulfil. And finally, I want to thank to my engaged Cristina since without her constant support and patience listening to me talking about “triangulitos”, I would not have been able to finish this thesis. Resumen La mayor´ıa de los problemas en f´ısica e ingenier´ıa se modelan mediante ecuaciones en derivadas parciales. Algunos de los modelos m´as representativos ser´ıan, por ejemplo, las ecuaciones de Maxwell, el problema de la elasticidad o las ecuaciones de Navier-Stokes. Generalmente, la resoluci´on de las ecuaciones en derivadas parciales se basa en el uso de m´etodos num´ericos que comienzan con un proceso de discretizaci´on y finalizan con la resoluci´on de un gran sistema de ecuaciones algebraicas que puede llegar a ser muy costosa. El desarrollo de m´etodos eficientes para esta ´ultima tarea es muy importante pues repercute en la posibilidad de realizar m´as ensayos con un menor costo computacional, y es el principal objetivo de esta tesis. Hay dos grandes grupos de m´etodos num´ericos para la resoluci´on de estos sistemas de ecuaciones, los m´etodos basados en los subespacios de Krylov, y los m´etodos multimalla. En aplicaciones especialmente dif´ıciles, es habitual combinar ambas metodolog´ıas, ya sea usando un m´etodo multimalla como precondicionador de un m´etodo basado en un subespacio de Krylov o este ´ultimo como suavizador en el m´etodo multimalla. Esta tesis se va a centrar en los m´etodos multimalla. Los m´etodos multimalla (4), que son conocidos por ser m´etodos muy r´apidos para la resoluci´on de problemas de tipo parab´olico y el´ıptico (49), se pueden clasificar en dos familias. Los geom´etricos que tienen un bajo consumo de memoria, pero solo pueden aplicarse sobre dominios regulares y mallas estructuradas, y los m´etodos multimalla algebraicos, que requieren m´as memoria pero se pueden aplicar a dominios de cualquier complejidad. Una alternativa que pretende obtener lo mejor de ambos m´etodos es el uso de m´etodos multimalla geom´etricos en mallas semi-estructuradas. Esta se basa en considerar una malla inicial no estructurada que se adapte a la geometr´ıa del dominio, para posteriormente refinar regularmente esta malla inicial hasta obtener una malla suficientemente fina en la que podamos obtener una soluci´on con la precisi´on deseada. iii De esta manera, obtenemos una malla que se adapta al dominio, y que es localmente estructurada, donde podemos aplicar m´etodos multimalla geom´etricos. Se han aplicado estos m´etodos sobre mallas rectangulares y triangulares, en estas ´ultimas considerando discretizaciones centradas en los v´ertices (20, 21). Las mallas de Delaunay son ampliamente utilizadas para la triangulaci´on de dominios debido a sus buenas propiedades. La malla dual de una triangulaci´on de Delaunay es conocida como malla de Voronoi. Entre otras buenas propiedades, se cumple que las rectas que unen los nodos de la malla de Voronoi son perpendiculares a los lados de la triangulaci´on de Delaunay. Por tanto, su uso para la discretizaci´on de ecuaciones en derivadas parciales es muy interesante, ya que para discretizar las derivadas en la direcci´on normal basta con un cociente de diferencias. Vamos a considerar triangulaciones acut´angulas, por lo que el punto de Voronoi siempre estar´a dentro del tri´angulo correspondiente. Sin embargo, su posici´on cambia en funci´on de la forma del tri´angulo, lo cual puede dar lugar a grandes anisotrop´ıas debidas al mallado y no al problema. Uno de los componentes m´as importantes de los m´etodos multimalla es el suavizador. El estudio de suavizadores para discretizaciones centradas en celdas es escaso, y nulo para mallas semi-estructuradas. Por tanto, el dise˜no de nuevos suavizadores es uno de nuestros principales objetivos. Para ello, deberemos estudiar las caracter´ısticas del problema y dise˜nar nuevos suavizadores que sean capaces de tener en cuenta dichas caracter´ısticas as´ı como de tratar con las posibles anisotrop´ıas inducidas por la malla. A la vez que desarrollamos un m´etodo multimalla gen´erico para dichas discretizaciones, nos centraremos en la resoluci´on del modelo Darcy-Oberbeck-Boussinesq (42), el cual se puede escribir de forma adimensional como sigue: ∇·u= 0,(1) u=−(∇P−Cz),(2) ∂C ∂t =−u·∇C+1 Ra∆C, (3) donde zes el vector (0,−1) en coordenadas cartesianas, y Ra denota el par´ametro de Rayleigh, el cual establece la dominancia de la difusi´on o la convecci´on. Comenzaremos con el dise˜no de un m´etodo multimalla para discretizaciones de problemas escalares en mallas de Voronoi asociadas a triangulaciones semi-estructuradas, m´as concretamente consideraremos el problema de Laplace. Tras obtener un multimalla eficiente para este problema, iremos a˜nadiendo complejidad al modelo, materiales no iv List of Figures 2.1 Unstructured mesh and its associated Voronoi grid. . . . . . . . . . . . . 7 2.2 Notation for neighbouring Voronoi points on an unstructured grid. . . . 8 2.3 (a) New basis in R2fitting the geometry of a uniform triangular grid, and local numeration for the regular Delaunay grid obtained on a triangular domain. (b) Corresponding Voronoi mesh. . . . . . . . . . . . . . . . . . 9 2.4 Stencils corresponding to two different grid-points: one associated with an up-oriented triangle and the other with a down-oriented triangle. . . 10 2.5 Voronoi mesh split into two sub-grids, and corresponding values of δi and δ0 i...................................... 11 2.6 Notation for neighbouring Voronoi points on a structured grid, characterized by angles αand β. .......................... 12 2.7 Variability of the ratios between the coefficients of a stencil associated with triangles with (a) α= 60oor (b) α= 88o,and βvarying among all its possible values satisfying the requirement of an acute triangulation. 14 2.8 Nested (left) and non-nested (right) hierarchies. . . . . . . . . . . . . . . 16 2.9 Restriction operator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.10 (a) Green-diamond smoother and (b) Red-wormy smoother. . . . . . . . 19 2.11 (a) Comparison of smoothers on an equilateral triangular grid. (b) History of the convergence for different numbers of refinement levels by using ω−red-blacksmoother. ........................ 21 2.12 (a) Comparison of smoothers on an almost-right triangular grid. (b) History of the convergence for different numbers of refinement levels by using diamond smoother. . . . . . . . . . . . . . . . . . . . . . . . . . . 22 xi LIST OF FIGURES 2.13 Experimentally computed convergence factors for the diamond smoother based multigrid and four smoothing steps, for different triangles in function of two of their angles. . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.14 (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. . . . . . . . . . . . . . . . . 24 2.15 Guideline to choose suitable smoothers to reach an asymptotic convergence factor about 0.1 on different triangles. . . . . . . . . . . . . . . . . 25 2.16 The Local Fourier Analysis and the multigrid code should provide the same asymptotic convergence rate. . . . . . . . . . . . . . . . . . . . . . 26 2.17 Infinite triangular Delaunay grid with the associated infinite Voronoi mesh and the basis fitting the grid . . . . . . . . . . . . . . . . . . . . . 27 2.18 Update order in Gauss-Seidel smoother for a Voronoi point associated with a down-oriented triangle (left) and a Voronoi point associated with an up-oriented triangle (right). . . . . . . . . . . . . . . . . . . . . . . . 31 2.19 Voronoi anisotropy appearing when almost-right triangular grids are considered. .................................... 32 2.20 (a) Diamond smoother and (b) wormy smoother on infinite grids. . . . . 32 2.21 Experimentally computed asymptotic convergence factors of a V-cycle by using different numbers of levels, for three different triangular configurations..................................... 36 2.22 (a) Coarsest unstructured grid. (b) Grid obtained after one regular refinementlevel.................................. 37 2.23 Different smoothers for the triangles composing the initial triangulation oftheA-shapeddomain............................ 38 2.24 Coarsest unstructured grid together with the associated Voronoi mesh. . 40 2.25 Different smoothers considered on each triangular block of the input grid. 41 2.26 A triangle structure stored in an array, keeping the row and column coherence ................................... 43 2.27 A triangle structure stored in two arrays, one for up-oriented triangles and another for down-oriented triangles. . . . . . . . . . . . . . . . . . . 44 2.28 A triangle structure stored in two arrays, one for up-oriented triangles and another for down-oriented triangles, with boundary data storage, B. 45 xii LIST OF FIGURES 2.29 (a) Initial unstructured grid of a domain. (b) Semi-structured grid obtained after one refinement level. . . . . . . . . . . . . . . . . . . . . . . 46 2.30 Three isolated triangles after a virtual splitting. . . . . . . . . . . . . . . 46 2.31 Overlapping update process. . . . . . . . . . . . . . . . . . . . . . . . . . 47 2.32 (a) Red-Black smoother in two triangles keeping the coherence between them. (b) Practical red-black smoother. . . . . . . . . . . . . . . . . . . 47 2.33 Error after two smoothing steps, following the relaxation strategy presentedinFigure2.15. ............................ 48 2.34 Wormy smoother extended to two lines of two different triangles at the sametime.................................... 49 2.35 Error after two smoothing steps, following the relaxation strategy presented in Figure 2.15, and applying the extra-relaxation process. . . . . 50 2.36 Extended wormy smoother, with an overlap of one node. . . . . . . . . . 51 2.37 Four nodes deep extra-relaxation. . . . . . . . . . . . . . . . . . . . . . . 52 3.1 Notation used to construct the stencil on a Voronoi point at (a) an uporiented triangle or at (b) a down-oriented triangle. . . . . . . . . . . . . 57 3.2 (a) Coarsest unstructured mesh for the first test case, and distribution of diffusion coefficients: κ= 0.333 ×105at the yellow region and κ= 2 at the green part. (b) Different smoothers for the triangles of the coarsest grid: white corresponds to red-black smoother, diamond smoother is represented by red, and wormy smoother by blue. . . . . . . . . . . . . . 59 3.3 (a) Coarsest unstructured mesh for the second test case, and distribution of diffusion coefficients: κ= 0.333×105at the yellow region and κ= 2 at the green part. (b) Different smoothers for the triangles of the coarsest grid: white corresponds to red-black smoother, diamond smoother is represented by red, and wormy smoother by blue. . . . . . . . . . . . . . 60 3.4 (a) Coarsest unstructured mesh and distribution of diffusion coefficients for the second experiment. Yellow color represents κ= 0.001 and green κ= 1. (b) Different smoothers for the triangles of the coarsest grid: white corresponds to red-black, diamond smoother is represented by red and wormy smoother by blue. . . . . . . . . . . . . . . . . . . . . . . . . 61 xiii LIST OF FIGURES 3.5 Comparison between direct discretization and Galerkin approach for different numbers of refinement levels. . . . . . . . . . . . . . . . . . . . . . 62 3.6 Streamlines of the circulatory flow in the square domain. . . . . . . . . . 68 3.7 (a) Grid after two refinement levels (b) Different smoothers for the triangles of the coarsest grid: white corresponds to red-black, diamond smoother is represented by red and wormy smoother by blue. . . . . . . 73 4.1 Unstructured mesh and its associated Voronoi grid with dots for the scalar values, and squares denoting the places where the projection of the vector values are stored. . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.2 (a) Restriction for a L1 type node (circles), where the red nodes have a weight of 1 4and the green ones of 1 12. (b) Prolongation for a L1 type node, where some nodes are updated by injection, whereas others, like the orange node, need to be calculated from the three coarse grid surrounding nodes. ..................................... 78 4.3 Vanka red-black smoother for staggered triangular grids. . . . . . . . . . 79 4.4 Process of a distributive smoother. . . . . . . . . . . . . . . . . . . . . . 81 4.5 History of the convergence for different smoothers in an equilateral triangular grid composed of 65536 triangles. . . . . . . . . . . . . . . . . . 83 4.6 History of the convergence for different smoothers in an almost righttriangular grid composed of 65536 triangles. . . . . . . . . . . . . . . . . 84 4.7 History of the convergence for the wormy smoother with different number ofrefinementlevels. ............................. 84 4.8 History of the convergence for the distributive smoother for different isosceles triangular grids. . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.9 History of the convergence for the distributive smoother for different numbers of refinement levels. . . . . . . . . . . . . . . . . . . . . . . . . 85 4.10 Guideline to choose suitable smoothers to reach an asymptotic convergence factor about 0.1 for different triangular grids characterized by two angles, αand β. ............................... 86 4.11 Different smoothers for the triangles of the coarsest grid: diamond smoother is represented by red, and wormy smoother by blue. . . . . . . . . . . . 87 5.1 Flowchart showing the fixed point method. . . . . . . . . . . . . . . . . 97 xiv LIST OF FIGURES 5.2 Initial unstructured meshes. . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.3 CO2concentration after 0.53 dimensionless seconds for the two consideredmeshes................................... 101 5.4 CO2concentration after 1.02 dimensionless seconds for the two consideredmeshes................................... 101 5.5 CO2concentration after 1.52 dimensionless seconds for the two consideredmeshes................................... 101 5.6 CO2concentration after 2.01 dimensionless seconds for the two consideredmeshes................................... 102 5.7 CO2concentration after 2.51 dimensionless seconds for the two consideredmeshes................................... 102 5.8 CO2concentration after 3.01 dimensionless seconds for the two consideredmeshes................................... 102 5.9 CO2concentration after 9.97 dimensionless seconds for the two consideredmeshes................................... 103 5.10 (a) CO2concentration after 0.07 dimensionless seconds. (b) CO2concentration after 0.25 dimensionless seconds. . . . . . . . . . . . . . . . . 104 5.11 (a) CO2concentration after 0.39 dimensionless seconds. (b) CO2concentration after 0.68 dimensionless seconds. . . . . . . . . . . . . . . . . 104 5.12 (a) CO2concentration after 1.22 dimensionless seconds. (b) CO2concentration after 2.50 dimensionless seconds. . . . . . . . . . . . . . . . . 105 5.13 Different smoothers for the triangles of the coarsest grid: white corresponds to red-black, diamond smoother is represented by red and wormy smootherbyblue. .............................. 105 5.14 (a) CO2concentration after 0.03 dimensionless seconds. (b) CO2concentration after 1.00 dimensionless seconds. . . . . . . . . . . . . . . . . 106 5.15 (a) CO2concentration after 2.08 dimensionless seconds. (b) CO2concentration after 2.51 dimensionless seconds. . . . . . . . . . . . . . . . . 106 5.16 (a) CO2concentration after 5.02 dimensionless seconds. (b) CO2concentration after 9.99 dimensionless seconds. . . . . . . . . . . . . . . . . 107 5.17 (a) CO2concentration after 1.02 dimensionless seconds. (b) CO2concentration after 2.08 dimensionless seconds. . . . . . . . . . . . . . . . . 107 xv LIST OF FIGURES 5.18 (a) CO2concentration after 3.00 dimensionless seconds. (b) CO2concentration after 4.02 dimensionless seconds. . . . . . . . . . . . . . . . . 108 5.19 (a) CO2concentration after 6.07 dimensionless seconds. (b) CO2concentration after 8.02 dimensionless seconds. . . . . . . . . . . . . . . . . 108 5.20 (a) Rayleigh number values in the considered domain. Blue stands for a Rayleigh number of 1000 and red for a number of 8000. (b) Different smoothers for the scalar case: white corresponds to red-black, diamond smoother is represented by red and wormy smoother by blue. . . . . . . 109 5.21 (a) CO2concentration after 1.09 dimensionless seconds. (b) CO2concentration after 2.47 dimensionless seconds. . . . . . . . . . . . . . . . . 109 5.22 (a) CO2concentration after 3.56 dimensionless seconds. (b) CO2concentration after 5.01 dimensionless seconds. . . . . . . . . . . . . . . . . 110 5.23 (a) CO2concentration after 7.02 dimensionless seconds. (b) CO2concentration after 9.99 dimensionless seconds. . . . . . . . . . . . . . . . . 110 5.24 (a) Rayleigh numbers in the considered domain. Light green represents a Rayleigh number of 100, Blue stands for 1000, red for 5000 and yellow for 8000. (b) Different smoothers for the scalar case: white corresponds to red-black, diamond smoother is represented by red and wormy smoother byblue. .................................... 111 5.25 (a) CO2concentration after 0.33 dimensionless seconds. (b) CO2concentration after 0.83 dimensionless seconds. . . . . . . . . . . . . . . . . 112 5.26 (a) CO2concentration after 1.17 dimensionless seconds. (b) CO2concentration after 2.33 dimensionless seconds. . . . . . . . . . . . . . . . . 112 5.27 (a) CO2concentration after 3.33 dimensionless seconds. (b) CO2concentration after 4.99 dimensionless seconds. . . . . . . . . . . . . . . . . 112 5.28 (a) CO2concentration after 6.66 dimensionless seconds. (b) CO2concentration after 9.99 dimensionless seconds. . . . . . . . . . . . . . . . . 113 A.1 Gauss-Seidel smoothing effect on the error. . . . . . . . . . . . . . . . . 124 A.2 Typesofcycles................................. 126 A.3 High and low frequencies, grey and white areas respectively, for standard coarsening. .................................. 129 xvi List of Tables 2.1 Number of cycles necessary to reduce the initial residual in a factor of 10−10,by using an F-cycle and the corresponding CPU-times. . . . . . . 24 2.2 LFA two-grid convergence factors, ρ2g, and measured two-grid convergence rates ρhfor an equilateral triangular grid. . . . . . . . . . . . . . . 33 2.3 LFA two-grid convergence factors, ρ2g, and measured two-grid convergence rates ρhfor an almost-right triangular grid. . . . . . . . . . . . . . 34 2.4 LFA two-grid convergence factors, ρ2g, and measured two-grid convergence rates ρhfor an almost-right triangular grid. . . . . . . . . . . . . . 35 2.5 LFA three-grid convergence factors, ρ3g, and experimentally computed factors using three and nine multigrid levels (the number of nodes is the same)...................................... 36 2.6 Number of iterations to reduce the initial residual in a factor of 10−10, corresponding asymptotic convergence rates and CPU-times for different numbers of refinement levels, by using a V-cycle and an F-cycle. . . . . 39 2.7 Number of iterations necessary to reduce the initial residual in a factor of 10−10,corresponding asymptotic convergence rates and CPU-times for different numbers of refinement levels, by using a V-cycle and an F-cycle. 42 2.8 Convergence rates using the modified wormy and a the standard one, with different relaxation parameters and for different very sharp isosceles triangles. ................................... 51 2.9 Convergence rates using the special extra-relaxation method for three very sharp isosceles triangles. . . . . . . . . . . . . . . . . . . . . . . . . 52 xvii LIST OF TABLES 3.1 Number of iterations to reduce the initial residual in a factor of 10−10 for both test cases: the inner rhombus subdomain (Figure 3.2) and the inner circle subdomain (Figure 3.3). . . . . . . . . . . . . . . . . . . . . 61 3.2 Number of iterations to reduce the initial residual in a factor of 10−10, corresponding asymptotic convergence rates (between brackets) and CPUtime for different numbers of refinement levels, by using a V-cycle and an F-cycle, and different numbers of smoothing steps ν= 2,3,4.The symbol −indicates divergence of the method. . . . . . . . . . . . . . . . 63 3.3 Convergence factors for different smoothers using a V(1,1)-cycle with different flow directions. . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.4 Number of iterations to reduce the initial residual in a factor of 10−10, corresponding asymptotic convergence rates and CPU-times for different numbers of refinement levels, using a V-cycle and an F-cycle. . . . . . . 68 3.5 Number of iterations to reduce the initial residual in a factor of 10−10, corresponding asymptotic convergence rates and CPU-times for different numbers of refinement levels, and by using only wormy smoothers. . . . 69 3.6 Asymptotic convergence rates for an almost right triangular grid characterized by angles (45, 85), and for an isosceles triangular grid characterized by angles (85,85), by using four different smoothers. The time discretization is τ= 0.01 and a grid with 16384 nodes is considered. . . 71 3.7 Convergence factors for different smoothers using a V(1,1)-cycle with different flow directions and τ= 0.01. ................... 72 3.8 Number of iterations to reduce the initial residual in a factor of 10−10, corresponding asymptotic convergence rates and CPU-times for different numbers of refinement levels, by using V-cycle and F-cycle, where τ= 0.01.73 3.9 Number of iterations to reduce the initial residual in a factor of 10−10, corresponding asymptotic convergence rates and CPU-times for different numbers of refinement levels, by using V-cycle and F-cycle, where τ= 0.00001..................................... 74 4.1 CPU time to perform 100 iterations of the relaxation scheme by using a grid composed of 65536 triangles. . . . . . . . . . . . . . . . . . . . . . . 82 xviii LIST OF TABLES 4.2 Number of cycles, asymptotic convergence factors and CPU times for different type of cycles and different numbers of refinement levels, by using the global distributive smoother. . . . . . . . . . . . . . . . . . . . 88 4.3 Number of cycles, asymptotic convergence factors and CPU times for different type of cycles and different numbers of refinement levels, by using the smoothers depicted in Figure 4.11. . . . . . . . . . . . . . . . 88 xix LIST OF TABLES xx 2.1 Discretization of an homogeneous diffusion problem on Voronoi grids fall outside their corresponding triangle. A Voronoi mesh asociated with a triangulation T, see Figure 2.1, is defined by the cenFigure 2.1: Unstructured mesh and its associated Voronoi grid. tres of the circumscribed circles of the triangles in T. Moreover, the segment connecting two neighboring Voronoi points goes trough the midpoint of the common edge of their corresponding triangles, being also perpendicular to it. This latter property allows us to approximate the derivative in the normal direction, by the difference quotient on both Voronoi points. Using the divergence theorem on a triangle Tof the triangulation T,the following balance equation holds: −Zl1∇v·n1dl1−Zl2∇v·n2dl2−Zl3∇v·n3dl3=ZT f(x) dx.(2.3) where niis the unit outward normal vector to the corresponding edge liof triangle T. Now, each of these line integrals is approximated as the length of the corresponding edge multiplied by the flux evaluated in the midpoint of the edge. Afterwards, we approximate such fluxes using the Voronoi points. With regard to the integral in the right-hand side, we consider the following approximation: ZT f(x) dx≈meas(T)f(xc),(2.4) where xcdenotes the Voronoi point of triangle T, and where meas(T) is the area of T. Denoting x1,x2,x3the Voronoi points of the triangles adjacent to Twith common edges l1, l2and l3,respectively, and dithe distance between points xcand xi, 7 2. SCALAR MULTIGRID METHODS ON SEMI-STRUCTURED VORONOI GRIDS Figure 2.2: Notation for neighbouring Voronoi points on an unstructured grid. di= dist(xc,xi),with i= 1,2,3,(see Figure 2.2), we finally obtain the equation corresponding to node xc: −1 meas(T) 3 X i=1 meas(li)vh(xi)−vh(xc) di=f(xc),(2.5) where vhis the grid function approximating the solution of the continuous problem on the Voronoi points. This scheme is first order accurate in a discrete H1−norm, whereas it is computationally observed that second order is achieved in L2−norm. The theoretical analysis of this approach can be seen in (26, 41, 52). 2.1.2 Discretization on structured triangular grids Now, we are going to consider the particular case of the discretization of problem (2.1) on a structured triangular grid. In this kind of grids, it is very usual to work in stencil notation because it takes advantage of the structured ordering of the unknowns which contribute in the discretization of a fixed grid-point. In a structured grid, any point is surrounded by the same grid-pattern, and using a suitable numbering of the gridpoints it is easy to capture this pattern in a small matrix or “stencil” which stores the contributions of the neighbouring unknowns. Then, first of all, a suitable numbering of the grid-points is needed. In triangular grids, a unitary basis of R2,{e1,e2},where e1, and e2are unit vectors defining the oblique coordinate system, is considered fitting the geometry of the triangle, as can be seen in Figure 2.3 (a). Hence, a local numeration 8 2.1 Discretization of an homogeneous diffusion problem on Voronoi grids (a) (b) Figure 2.3: (a) New basis in R2fitting the geometry of a uniform triangular grid, and local numeration for the regular Delaunay grid obtained on a triangular domain. (b) Corresponding Voronoi mesh. can be fixed according to the definition of the spatial basis. In this way, a manner of numbering nodes very convenient for identifying the neighbouring nodes can be defined. We consider a triangular grid arising on a triangular domain by applying a fixed number of regular refinement steps `. This is done in the way that, on each refinement step every triangle is divided into four congruent ones by connecting the midpoints of their edges. Then, we can define the corresponding grid in the following way: G`={x=k1h1e1+k2h2e2|k1= 0,...,2`, k2= 0, . . . , k1},(2.6) where h= (h1, h2) is the grid spacing associated with the refinement level `(h1is the grid spacing in the direction of e1,and h2in the direction of e2), so that the grid G`can also be denoted by Gh. Thus, for a refinement level `, a local numeration with double index (k1, k2), k1= 0,...,2`,k2= 0, . . . , k1, is used in such a way that the indexes of the vertices of the triangle are (0,0), (2`,0), (2`,2`), as it can also be observed in Figure 2.3(a) for `= 2. On the other hand, the considered discretization is based on the dual Voronoi mesh, represented in Figure 2.3(b). In the particular case in which a structured grid as considered here is used, the obtained finite difference scheme becomes different depending on the grid-point. More concretely, one-half of the grid-points, those corresponding to an up-oriented triangle, have the same equation and the other half, those corresponding to a down-oriented triangle, have a “mirror image stencil”, see Figure 2.4. In this 9 2. SCALAR MULTIGRID METHODS ON SEMI-STRUCTURED VORONOI GRIDS Figure 2.4: Stencils corresponding to two different grid-points: one associated with an up-oriented triangle and the other with a down-oriented triangle. sense, the Voronoi mesh, denoted by Vhcould be split up into two sub-grids Vu h(associated with the up-oriented triangles) and Vd h,(corresponding to the down-oriented triangles), as seen in Figure 2.3(b). These sub-grids can be defined from the grid Gh, in the following way: Vu h={xu k1,k2= (k1+δ1)h1e1+ (k2+δ2)h2e2|k1h1e1+k2h2e2∈Gh},(2.7) Vd h={xd k1,k2= (k1+δ0 1)h1e1+ (k2+δ0 2)h2e2|k1h1e1+k2h2e2∈Gh},(2.8) where δi, δ0 i,with i= 1,2,are suitable scalar values to reach Voronoi points from the primal ones following the considered local coordinate system, see Figure 2.5. Then, a grid-function, vh,defined on the Voronoi mesh Vh,could be split up into two different sub-grid functions, vu hand vd h,associated with sub-grids Vu hand Vd h,respectively. In this way, given an arbitrary pair (k1, k2) associated with a node of Gh,the equations corresponding to the two Voronoi points xu k1,k2and xd k1,k2,are given by Luu hvu h(xu k1,k2) + Lud hvd h(xd k1,k2) = fu h(xu k1,k2),(2.9) Ldu hvu h(xu k1,k2) + Ldd hvd h(xd k1,k2) = fd h(xd k1,k2),(2.10) where these “scalar” operators are given in stencil form as: 10 2.1 Discretization of an homogeneous diffusion problem on Voronoi grids Figure 2.5: Voronoi mesh split into two sub-grids, and corresponding values of δiand δ0 i. Luu h=1 meas(T)       0 0 0 0 3 X i=1 li di 0 0 0 0       , Lud h=1 meas(T)       0 0 0 0−l1 d1−l3 d3 0−l2 d2 0       , Ldu h=1 meas(T)       0−l2 d2 0 −l3 d3−l1 d1 0 0 0 0       , Ldd h=1 meas(T)       0 0 0 0 3 X i=1 li di 0 0 0 0       . 2.1.2.1 Stencil depending on two angles characterizing the triangular grid. An important feature to take into account in a stencil is the strength of the connections between the involved unknowns. Each entry of the stencil defines the coefficient of the corresponding unknown in the equation associated with the central point. Then, if some coefficient is large relative to the other coefficients in the same equation, then a small change in the value of the associated unknown has more effect on the value of the unknown at the central point, and therefore, we will say that these both unknowns are strongly connected, see (7). In order to know a priori the strong and weak connections between neighbouring unknowns depending on the grid geometry, we are going to rewrite the stencils as a function of some parameters characterizing the grid, that is, two angles, αand β, and the length l, of one edge of an arbitrary triangle of the grid, see Figure 2.6. As we will 11 2. SCALAR MULTIGRID METHODS ON SEMI-STRUCTURED VORONOI GRIDS see, this is going to be very useful for the design of smoothers for different geometries, taking into account the strong connections appearing in the stencils. Therefore, we are going to describe in detail the computation of the stencil for a Voronoi grid-point, associated with an arbitrary down-oriented triangle, xd k1,k2in Vd h.Since all the parameters involved in the stencil (area, distances between Voronoi points and lengths of the edges) are independent of the chosen coordinate system, for simplicity, the coordinates of the points involved in such stencil can be computed in the Cartesian coordinate system with respect to the origin, see Figure 2.6. And, in terms of the previously explained geometric parameters, they result in the following: xd k1,k2=l 22 cos αsin β sin(α+β)+ 1,2 sin αsin β sin(α+β)+ cot(α+β),(2.11) xu k1,k2=l 2(3,−cot(α+β)) ,(2.12) xu k1−1,k2=l 2(1,−cot(α+β)) ,(2.13) xu k1,k2+1 =l 22 cos αsin β sin(α+β)+ 1,2 sin αsin β sin(α+β)−cot(α+β).(2.14) Figure 2.6: Notation for neighbouring Voronoi points on a structured grid, characterized by angles αand β. Due to the fact that the area of an arbitrary triangle Tis given in terms of the 12 2.1 Discretization of an homogeneous diffusion problem on Voronoi grids geometric parameters as meas(T) = l2sin αsin β 2 sin(α+β),(2.15) and the lengths of the sides of Tare l2=l, l1=lsin β sin(α+β),and l3=lsin α sin(α+β),after computing dii= 1,2,3 from (2.11)-(2.14), we can finally obtain the stencils: Ldu h=2 sin(α+β) l2sin αsin β  0 tan(α+β) 0 −tan α−tan β0 0 0 0  ,(2.16) Ldd h=2 sin(α+β) l2sin αsin β  000 0−tan(α+β) + tan α+ tan β0 000 .(2.17) As previously commented, for a Voronoi grid-point, associated with an arbitrary uporiented triangle, xu k1,k2in Vu hthe stencil would be the “mirror image stencil” of (2.16)- (2.17), that is, Luu h=2 sin(α+β) l2sin αsin β  000 0−tan(α+β) + tan α+ tan β0 000 ,(2.18) Lud h=2 sin(α+β) l2sin αsin β  0 0 0 0−tan β−tan α 0 tan(α+β) 0  .(2.19) Notice that depending on the angles characterizing the grid, some strong connections appear between some unknowns. In particular, in Figure 2.7 we show the variability of the stencil weights, for a Voronoi point xd k1,k2associated with a down-oriented triangle, depending on one of the angles of the triangulation. More concretely, we display the ratios between the three extra-diagonal coefficients in the equation and the central term, that is, sdu 00 /sdd 00, sdu 01 /sdd 00,and sdu −10/sdd 00,(see notation in Figure 2.4) for two different situations in which angle αis fixed and angle βvaries among all its possible values satisfying the requirement of an acute triangulation. In Figure 2.7(a), αis 60o, and we can see that when β= 60o,that is for an equilateral triangle, the strength of the three connections in the stencil is the same. When angle βis small or large enough, we observe that two of the connections are weak whereas the other one is very strong. On the other hand, in Figure 2.7(b), αis fixed as 88o.We can see that for almost the whole range of values of βone of the connections in the stencil is strong and the other 13 2. SCALAR MULTIGRID METHODS ON SEMI-STRUCTURED VORONOI GRIDS (a) (b) Figure 2.7: Variability of the ratios between the coefficients of a stencil associated with triangles with (a) α= 60oor (b) α= 88o,and βvarying among all its possible values satisfying the requirement of an acute triangulation. ones are much weaker, whereas when βis either very small or very large, then two of the connections in the stencil are of equal strength. These strong connections between the unknowns involved in a stencil will have to be taken into account in the design of the smoothers in a geometric multigrid method, as will be discussed in next section. 2.2 Multigrid method In this section we are going to present the multigrid components that we are going to use. For a basic explanation of multigrid methods, we refer the reader to Appendix A. 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 the following subsections the proposed cell-centred 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-centred 14 2.2 Multigrid method 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 semistructured 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 block-wise 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 chapter. 2.2.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. 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 fine scale to obtain the desired approximation of the solution. Once the hierarchy of grids is defined, we have to choose suitable inter-grid transfer operators as well as the corresponding discrete operators on the coarse grids. Regarding the coarse-grid approximation, a direct PDE discretization has been used. When vertex-centred discretizations are considered on triangular grids, grid-points lying on coarser grids also belong to the finer grids, giving rise to a so-called nested hierarchy of grids. However, when the considered cell-centred 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 2.8. 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 15 2. SCALAR MULTIGRID METHODS ON SEMI-STRUCTURED VORONOI GRIDS 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. Figure 2.8: Nested (left) and non-nested (right) hierarchies. 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 =   000 0 1/4 0 000   ,(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 =   000 0 1/4 0 000   , (2.20) as shown in Figure 2.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 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 coarsegrid correction part of the algorithm. Actually, the design of suitable smoothers for cell-centred grids is a challenge in this context. 2.2.2 Smoothers The smoother usually plays an important role in multigrid algorithms, mainly in the geometric approach. The choice of a suitable smoother is an important feature to 16 2.2 Multigrid method eleven cycles. Besides, in Figure 2.12(b), 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 behaviour can be seen in Figure 2.13, where asymptotic convergence factors of the diamond smoother based multigrid are shown for a wide range of pairs of angles (α, β) characterizing the grid. Figure 2.13: 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 block-wise 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 2.14(a), 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 2.14(b). 23 2. SCALAR MULTIGRID METHODS ON SEMI-STRUCTURED VORONOI GRIDS Equilateral Right-Triangle Sharp triangle (ω−red-black) (Diamond smoother) (Wormy smoother) Levels Unknowns Cycles Time (s) Cycles Time (s) Cycles Time (s) 6 4096 7 0.03 8 0.03 6 0.04 7 16384 7 0.09 9 0.10 6 0.13 8 65536 7 0.35 9 0.39 7 0.58 9 262144 7 1.35 9 1.55 7 2.28 10 1048576 7 5.76 9 6.81 7 9.66 Table 2.1: Number of cycles necessary to reduce the initial residual in a factor of 10−10, by using an F-cycle and the corresponding CPU-times. (a) (b) Figure 2.14: (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. From the results presented in this section, 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-right triangles, and finally the appropriate block collective wormy smoother when triangulations with a small angle appear. To summarize, in Table 2.1 we show for the three considered representative triangles the results corresponding to the best smoother for each geometry. In particular, for different numbers of refinement levels, the number of cycles necessary to reduce the initial residual in a factor of 10−10,and the CPU-time are shown. As we can observe, a convergence independent on the number of unknowns is obtained in all the three cases. 24 2.3 Local Fourier analysis From the practical point of view, for any given triangular geometry it would be nice to be capable of choosing a suitable smoother in order to reach a desired convergence factor. Moreover, for semi-structured grids, it is imperative to know the smoother to use for each triangle of the input grid to achieve a local desired convergence factor (45). In order to reach this, a set-up phase has been implemented in the multigrid code; it consists of reading an already calculated database containing the most efficient strategy depending on the angles of the triangle. That is, taking into account that wormy smoother is about twice expensive than ω−red-black smoother, and diamond smoother is about ten per cent more expensive than this latter, we choose for each triangle the cheapest smoother which gives a convergence factor below a desired one. Following this strategy, the corresponding guideline to reach a global convergence factor about 0.1 is shown in Figure 2.15. Figure 2.15: Guideline to choose suitable smoothers to reach an asymptotic convergence factor about 0.1 on different triangles. 2.3 Local Fourier analysis As we have seen in previous sections, the design of an efficient multigrid method may not be easy. The necessity of a good smoothness, and its interplay with the coarse-grid 25 2. SCALAR MULTIGRID METHODS ON SEMI-STRUCTURED VORONOI GRIDS correction is essential to achieve a satisfactory convergence rate. Up to this point, we have been working with numerical experiments, which means that we had to code the whole multigrid method in order to know whether it works or not. Taking into account the amount of work that this requires, it seems logical to try to know a priori if what we are going to perform is a good strategy. For that reason, the Local Fourier analysis, LFA, was initially developed by Brandt in (4) and afterwards extended by him in (6). It was focused on rectangular grids and recently, it was extended to triangular meshes in (20). This tool calculates the asymptotic convergence rate of a multigrid algorithm. To do that, it requires some assumptions: it neglects the effect of boundary conditions, by considering the discrete operator to be defined on an infinite grid, and also if the equation is non-linear or it has non-constant coefficients it must be locally linearized (49). Despite these limitations, this method is a very useful tool to design efficient multigrid methods for new problems. Also, it can be used to test if a multigrid code has any bugs, since LFA predicts the asymptotic convergence rate that the code should provide, see Figure 2.16. In Appendix A, we have included a basic introduction to the smoothing analysis. Figure 2.16: The Local Fourier Analysis and the multigrid code should provide the same asymptotic convergence rate. Local Fourier analysis cannot be straightforwardly applied to discretizations on regular Voronoi meshes, since two different stencils appear for different grid-points, depending on the orientation of the corresponding triangles. This latter, together with the fact that triangular meshes are considered, make this analysis different from the usual one. Regarding the extension to triangular meshes, the key is to write the Fourier transform using coordinates in non-orthogonal bases fitting the structure of the grid, as in Section 2.1.2. For local Fourier analysis, a regular infinite grid is assumed. With this purpose, we extend the definition of grid Ghin (2.6) to the corresponding infinite 26 2.3 Local Fourier analysis grid in the following way: Gh={xk1,k2=k1h1e1+k2h2e2|k1, k2∈Z},(2.25) as seen in Figure 2.17. In this way the associated Voronoi mesh is also extended to an infinite grid Vh,which can be decomposed again into two infinite sub-grids as Vh=Vu hSVd h.Now, we must extend the discrete problem to the whole infinite grid Vh. Figure 2.17: Infinite triangular Delaunay grid with the associated infinite Voronoi mesh and the basis fitting the grid However, for this purpose, we have to take into account that equations at points on Vu h and Vd hare different. Thus, we can define the application of the discrete operator to a grid-function vhon Vhas: Lhvh(x) =          X (i,j)∈Iuu suu i,j vh(xu k1+i,k2+j) + X (i,j)∈Iud sud i,jvh(xd k1+i,k2+j),x=xu k1,k2∈Vu h, X (i,j)∈Idu sdu i,jvh(xu k1+i,k2+j) + X (i,j)∈Idd sdd i,jvh(xd k1+i,k2+j),x=xd k1,k2∈Vd h, (2.26) where the coefficients in (2.26) are sdu 0,1=2 sin(α+β) l2sin αsin βtan(α+β), sdu −1,0=−2 sin(α+β) l2sin αsin βtan α, sdu 0,0=−2 sin(α+β) l2sin αsin βtan β, sdd 0,0=−(sdu 0,0+sdu 1,0+sdu 0,−1). (2.27) Due to the geometric relations between an upand a down-oriented triangle, it is easy to see that the coefficients satisfy sud 0,−1=sdu 0,1, sud 1,0=sdu −1,0, sud 0,0=sdu 0,0, suu 0,0=sdd 0,0. 27 2. SCALAR MULTIGRID METHODS ON SEMI-STRUCTURED VORONOI GRIDS and the subsets Iuu, Iud, Idu,and Idd give the connections of Voronoi points associated with up-oriented triangles with themselves and with those associated with downoriented triangles, and vice versa. In particular, in our case we have Iuu ={(0,0)}, Iud ={(0,0),(1,0),(0,−1)}, Idu ={(0,0),(−1,0),(0,1)},and Idd ={(0,0)}. A local Fourier analysis is based on the Fourier modes, which have the following form: ϕh(θ,x) = eıθ·x/h=eı(θ1x1/h1+θ2x2/h2),with θ= (θ1, θ2)∈Θh= [−π, π)2,and where x is a grid-point. For node-based problems, these Fourier modes turn out to be eigenvectors of grid operators that can be represented by a single stencil, see (4, 57). However, for the analysis that we want to develop, ϕh(θ,x) is not directly useful, and we have to define the corresponding “Fourier modes”, taking into account that we consider a cell-centered discretization on Vh,and the fact that the discrete operator cannot be represented by a single stencil. Then, we define the following grid-functions: φh(θ,x) = αuφu h(θ,x) + αdφd h(θ,x),x∈Vh, αu, αd∈C,(2.28) where φu h(θ,x) = ϕh(θ,x),x∈Vu h 0,x∈Vd h ,and φd h(θ,x) = 0,x∈Vu h ϕh(θ,x),x∈Vd h .(2.29) This splitting in (2.28) is a direct consequence of the fact that the stencils for the points on Vu hand Vd hare different. In this way, these grid-functions will play the same role in the LFA for the discretizations considered here as the Fourier modes ϕh(θ,x) in the LFA of standard node-based discretizations. Then, we define the following space of grid-functions: F(Vh) = {φh(θ,·) = αuφu h(θ,·) + αdφd h(θ,·), αu, αd∈C,θ∈Θh},(2.30) which will play the role of our “Fourier space”. In this way, we can prove that any discrete operator as in (2.26) leaves invariant subspace F(Vh).In fact, if we apply operator Lhto a grid-function φh∈F(Vh),depending on the grid-point and taking into account definitions in (2.7)-(2.8) the following holds: 28 2.3 Local Fourier analysis (Lhφh(θ,·))(xu k1,k2) = αuX (i,j)∈Iuu suu i,j φu h(xu k1+i,k2+j) + αdX (i,j)∈Iud sud i,jφd h(xd k1+i,k2+j) =eıθ·xu k1,k2/h αuX (i,j)∈Iuu suu i,j eı(θ1i+θ2j)+αdX (i,j)∈Iud sud i,jeı(θ1(i+δ0 1−δ1)+θ2(j+δ0 2−δ2)) , (Lhφh(θ,·))(xd k1,k2) = αuX (i,j)∈Idu sdu i,jφu h(xu k1+i,k2+j) + αdX (i,j)∈Idd sdd i,jφd h(xd k1+i,k2+j) =eıθ·xd k1,k2/h αuX (i,j)∈Idu sdu i,jeı(θ1(i+δ1−δ0 1)+θ2(j+δ2−δ0 2)) +αdX (i,j)∈Idd sdd i,jeı(θ1i+θ2j) . Then, from these expressions we obtain (Lhφh(θ,·)) = Lhφu hφd hαu αd = αuX (i,j)∈Iuu suu i,j eı(θ1i+θ2j)+αdX (i,j)∈Iud sud i,jeı(θ1(i+δ0 1−δ1)+θ2(j+δ0 2−δ2)) φu h + αuX (i,j)∈Idu sdu i,jeı(θ1(i+δ1−δ0 1)+θ2(j+δ2−δ0 2)) +αdX (i,j)∈Idd sdd i,jeı(θ1i+θ2j) φd h =φu hφd he Lh(θ)αu αd=φu hφd hβu βd,(2.31) from which the invariance property is demonstrated, and where e Lh(θ) =     X (i,j)∈Iuu suu i,j eı(θ1i+θ2j)X (i,j)∈Iud sud i,jeı(θ1(i+δ0 1−δ1)+θ2(j+δ0 2−δ2)) X (i,j)∈Idu sdu i,jeı(θ1(i+δ1−δ0 1)+θ2(j+δ2−δ0 2)) X (i,j)∈Idd sdd i,jeı(θ1i+θ2j)    , (2.32) is the representation of Lhin F(Vh).Notice that in this case, this representation is a 2×2−matrix. In the particular case of the discretization considered for our model problem, the Fourier domain representation of the discrete operator depending on the angles characterizing the triangular grid reads e Lh(θ) = 1 meas(T)tα+tβ−tαβ E(tαβe−ıθ2−tβ−tαeıθ1) E−1(tαβeıθ2−tβ−tαe−ıθ1)tα+tβ−tαβ (2.33) where the following abbreviations have been used: tα= tan(α), tβ= tan(β), tαβ = tan(α+β),and E=eı(θ1(δ0 1−δ1)+θ2(δ0 2−δ2)). 29 2. SCALAR MULTIGRID METHODS ON SEMI-STRUCTURED VORONOI GRIDS 2.3.1 Smoothers and their representations in the Fourier space In this work some classical smoothers will be considered together with some novel smoothers which have been appropriately designed to deal with the particular anisotropies that can appear in the Voronoi meshes. All these smoothers are based on a splitting of the discrete operator Lhas Lh=L+ h+L− h,that is, a decomposition on the positive and negative parts of the operator which correspond to the updated and non-updated unknowns before the current step, see (49). Next, each smoother is described and the corresponding Fourier domain representation is given. For the sake of simplicity, only the positive part of the operator is displayed. 2.3.1.1 Jacobi First, we consider the damped Jacobi iteration. In this case, L+ his just equal to the diagonal part of operator Lh,denoted by Dh,and the iteration matrix of this smoother is given by Sh=Ih−ωD−1 hLh,where Ihis the identity operator and ωis a relaxation parameter. The Fourier domain representation of operator Sh,e Sh(θ),is e Sh(θ) = e Ih(θ)−ωe Dh(θ)−1e Lh(θ),where e Ih(θ) is the 2 ×2−identity matrix, e Lh(θ) is given in (2.33), and e Dh(θ) = 1 meas(T)tα+tβ−tαβ 0 0tα+tβ−tαβ . 2.3.1.2 Gauss-Seidel Next, we consider the lexicographic Gauss-Seidel method, in which the grid-points are updated with an order dictated by the numbering of the Voronoi points. Since two Voronoi points are associated with the same pair of indexes (k1, k2),first the Voronoi point associated with a down-oriented triangle is relaxed and then that associated with the up-oriented triangle is updated. Figure 2.18 illustrates the update order. Considering this order, the Fourier representation of the positive part of the splitting is in this case e L+ h(θ) = 1 meas(T)tα+tβ−tαβ E(tαβe−ıθ2−tβ) −E−1tαe−ıθ1tα+tβ−tαβ .(2.34) 30 2.3 Local Fourier analysis Figure 2.18: Update order in Gauss-Seidel smoother for a Voronoi point associated with a down-oriented triangle (left) and a Voronoi point associated with an up-oriented triangle (right). 2.3.1.3 Red-black smoother Now, we consider a point-wise pattern relaxation that will be very efficient for equilateral triangulations. Since the Voronoi grid is naturally subdivided into two sub-grids Vu hand Vd h,we can consider a scheme which updates first the grid-points in Vu h,and secondly those points in Vd h.Thus, this scheme consists of two partial relaxation steps, in the first one the Voronoi points associated with up-oriented triangles are relaxed, and in the second one those grid-points associated with down-oriented triangles are updated. In this case, the Fourier representation of the positive part of the decomposition reads e L+ h(θ) = 1 meas(T)(tα+tβ−tαβ)/ωu0 E−1(tαβeıθ2−tβ−tαe−ıθ1) (tα+tβ−tαβ)/ωd,(2.35) where ωuand ωdare relaxation parameters associated with each one of the partial steps. 2.3.1.4 Diamond smoother Here we introduce a novel block-wise smoother suitable for one type of the Voronoi anisotropic grids arising when almost-right triangular grids are considered, see Figure 2.19. It consists of simultaneously relaxing two Voronoi points, each one associated with a different oriented triangle. This smoother will be efficient when these two points are very close. For instance, in the case of the almost-right triangulation displayed in Figure 2.19, the coupling of the Voronoi points will be as shown in Figure 2.20 (a), and then the Fourier domain representation of the positive part of the operator is 31 2. SCALAR MULTIGRID METHODS ON SEMI-STRUCTURED VORONOI GRIDS Figure 2.19: Voronoi anisotropy appearing when almost-right triangular grids are considered. e L+ h(θ) = 1 meas(T)tα+tβ−tαβ E(tαβe−ıθ2−tβ−tαeıθ1) −E−1tβtα+tβ−tαβ .(2.36) (a) (b) Figure 2.20: (a) Diamond smoother and (b) wormy smoother on infinite grids. 2.3.1.5 Wormy smoother We finish the description of the smoothers with the so-called wormy-smoother. As in the previous case, this is a block-wise smoother, and as will be seen in the results section it will be very efficient for the remaining case of Voronoi anisotropic meshes, that is, in the case of triangulations with a very small angle. This smoother consists of simultaneously updating all the Voronoi points associated with the triangles situated 32 2.4 Numerical experiments on semi-structured grids V(2,2) F(2,2) Levels Unknowns Cycles ρhTime (s) Cycles ρhTime (s) 4 51456 9 0.16 1.78 8 0.09 3.27 5 205824 9 0.19 3.62 8 0.11 6.32 6 823296 10 0.21 10.16 9 0.11 15.96 7 3293184 11 0.23 34.69 9 0.12 43.66 8 13172736 11 0.26 123.65 9 0.12 142.52 Table 2.6: Number of iterations to reduce the initial residual in a factor of 10−10,corresponding asymptotic convergence rates and CPU-times for different numbers of refinement levels, by using a V-cycle and an F-cycle. 2.4.2 Convection-diffusion problem on a square domain The strategy proposed can be applied to more complex problems as, for example convection-diffusion, which can be written as: −∇·(∇v+bv) = f, in Ω,(2.38) where b(x) is a given velocity field, whose divergence is assumed to be zero. In order to obtain a difference scheme by the cell-centred finite volume method, we follow the same approach that we have explained in detail in Section 2.1, by using a central difference scheme to approximate the convective term, see (41). In this numerical experiment an square domain of unit length and Dirichlet boundary conditions are considered, and a constant vector b= (1,0) is fixed in the whole domain. Thus, the following equation on each of the grid-nodes xcresults: −1 meas(T) 3 X i=1 meas(li)vh(xi)−vh(xc) di +b·ni vh(xi) + vh(xc) 2=f(xc). (2.39) We consider an initial unstructured grid, composed of 96 triangles, as seen in Figure 2.24, in which, for illustration, the dual Voronoi mesh has been displayed. The hierarchy of grids is obtained by regular refinement. As the convective part of the problem is not dominant, and its derivatives are of lower order, the behaviour of the multigrid will be similar to that obtained for a pure diffusive problem. Therefore we will follow the guideline given in Figure 2.15 to choose the suitable local smoother on each input triangle, and this selection is displayed in Figure 2.25. The proposed 39 2. SCALAR MULTIGRID METHODS ON SEMI-STRUCTURED VORONOI GRIDS Figure 2.24: Coarsest unstructured grid together with the associated Voronoi mesh. geometric multigrid method is applied to solve the corresponding large sparse linear system of equations. First, an F(2,2)−cycle is used to test the independence of the multigrid convergence with regard to the discretization parameters. In Table 2.7, for different numbers of refinement levels, the asymptotic convergence rate, ρh,and the number of iterations necessary to reduce the initial residual in a factor of 10−10,are displayed together with the CPU-time. Similarly to the previous numerical experiment, the results corresponding to apply a V(2,2)−cycle are also shown in this table. Again, a deterioration of the asymptotic convergence factor is observed, but the application of this cycle can be interesting due to the small number of iterations necessary to reach the convergence. 2.5 Implementation In this section, we are going to go deeper into the details of programming a multigrid algorithm for cell-centred discretizations on Voronoi meshes associated with semistructured triangular grids. The semi-structured approach requires the use of a connectivity array for the unstructured triangulation. Nevertheless, the triangles resulting from the regular refine40 2.5 Implementation Figure 2.25: Different smoothers considered on each triangular block of the input grid. ment of this initial unstructured mesh, have a structure ordering and therefore, the connectivity matrix is unnecessary. Hence, its use for all the triangles would result in a waste of resources, and stencils are used for that purpose. However, storing information about triangles is a challenge by itself since the data structures use to have a square shape. Also, we will comment how to overcome the problems that appear in the application of multigrid on semi-structured grids, by using extra-relaxation processes. And more specifically, the particular case in which the domain is composed of two extremely sharp triangles connected. 2.5.1 Data storage for cell-centred discretizations on structured triangular grids In all the numerical methods for the solution of partial differential equations, it is necessary to read data from the memory several times, apart from performing many mathematical operations. While the mathematical operations are done extremely fast, the access to the memory is slow and strongly depends on the data structure. Therefore, the way in which we store the information as well as the speed to read and write data from the memory will be a key point in the performance of the overall code. Thus, 41 2. SCALAR MULTIGRID METHODS ON SEMI-STRUCTURED VORONOI GRIDS V(2,2) F(2,2) Levels Unknowns Cycles ρhTime (s) Cycles ρhTime (s) 4 24576 9 0.19 0.74 8 0.09 1.19 5 98304 10 0.20 1.76 8 0.09 2.54 6 393216 10 0.21 4.54 8 0.09 6.07 7 1572864 11 0.22 15.53 8 0.09 17.11 8 6291456 11 0.22 55.16 8 0.09 56.33 Table 2.7: Number of iterations necessary to reduce the initial residual in a factor of 10−10,corresponding asymptotic convergence rates and CPU-times for different numbers of refinement levels, by using a V-cycle and an F-cycle. we can improve or reduce the performance of the algorithm depending on the used data structure. A desired data structure must be easy to use by the user and suitable for the performance of the computer. Therefore, one has to design a user friendly data structure depending also on the peculiarities of the programming language. In the particular case of this thesis, Fortran90 was the chosen language. In Fortran, the natural order of storage is column-major. Not following this order may make the data access time three times slower (11). Many programmers that need to solve PDEs claim that the most efficient way to store the data is by using a simple one-dimensional array. This is commonly used when rectangular grids are considered. In this case, each row has the same number of nodes, which makes easy the use of one dimensional arrays. To access to the data in column x, one just need to read from the memory position (x−1) ∗(size(column)) + 1 to x∗size(column).However, in the case of triangles, the size of the column is not fixed. Therefore, to obtain the information from a certain column, one needs to use a formula instead of a constant to know the beginning and the end of that column. Due to the already mentioned difficulties, we propose a simpler data structure based on 2-dimensional arrays. The storage of triangles have the complication that, if one wants to keep the logical structure of the triangle, this means that each row of a triangle is in a different column of the array, that array would be half empty, see Figure 2.26. This way of storage leaves a lot of empty space, and therefore it is undesirable. Nonetheless, if we add an extra row in the upper part of the array, then we have the same number of empty places, where we can store, for example, the right hand side of 42 2.5 Implementation Figure 2.26: A triangle structure stored in an array, keeping the row and column coherence the equation. However it is still not easy to move through nodes. For example, triangle number two is just below triangle number eight but in the array that coherence is not kept. Besides, to access to the right hand side data there is no a simple way. For example, node 1, that would be stored in the position (1,2) of the array, would have its right hand side in position (4,7).Despite possible, this way of storing is not easy and therefore likely to make mistakes. Hence, we will consider a new data structure. If we have a look to Figure 2.26, we can distinguish two sorts of triangles. Some of them are up-oriented triangles and the rest are down-oriented triangles. Each uporiented triangle is surrounded only by down-oriented triangles and vice-versa. This means, that to access to the neighbouring triangles the up-oriented triangles will only need information from the down-oriented triangles and the other way around. In this way, a new possible data structure, a modification of the previous one, arises in a natural way, to store each kind of triangles in a different array. However, we still have the problem of the inherent difficulty to access to the right hand side data. Regarding this latter, one solution is to keep the diagonal of the array empty, and to store the data of the right hand side in a symmetrical fashion, making simpler the access to the right hand side data, since the solution and the right-hand side data corresponding to the same node will have a symmetric numbering. For example: following Figure 2.27 the position of node 1 is (1,2) and its right hand side is stored in (2,1),which is much 43 2. SCALAR MULTIGRID METHODS ON SEMI-STRUCTURED VORONOI GRIDS easier to find than in the previous case, were both position had nothing in common. Also, the access to those triangles that are surrounding a concrete one is very easy. For example, if we want to find the neighbours of the up-oriented triangle number 10,whose storage position is (2,4),we just have to go to the array of down-oriented triangles and get the nodes corresponding to (2,4),(2,3) and (1,3), nodes 11,9 and 4 respectively. That is, the way to access to the neighbours is more or less the same as graphically we can observe, the left triangle is also in the left position in the array, the down triangle is also below and only in the right node one has to remember that it is stored in the same position of the other array. Similarly, for a down-oriented triangle, the method to obtain the neighbours is symmetric. That is, considering triangle 9 in the position (2,3) of the down-oriented triangles array, the neighbours, 8 (left), 10 (right) and 13 (up), are stored in the up-oriented triangles array in positions (2,3),(2,4) and (3,4) respectively. Figure 2.27: A triangle structure stored in two arrays, one for up-oriented triangles and another for down-oriented triangles. On the other hand, there is still memory that is not going to be used. More concretely, the amount of useless memory is exactly of the size of the length of the array, lenght(array), whereas the total amount of data is lenght(array)2. Hence, the memory that we are not going to use is almost irrelevant compared with the total amount of used memory. Regarding the storage of data corresponding to Dirichlet boundary conditions, the idea is to modify the down-oriented triangles array, by adding space to store the boundary data. We should store it in the places where we would look for neighbours from the up-oriented triangles that are touching the boundaries. For that, it is necessary to make the down-oriented array bigger, then it has to start from node (0,0),while the 44 2.5 Implementation up-oriented one must begin from (1,1).In this way, taking into account the boundary data will be as natural as accessing the neighbouring triangles. For example, for triangle 1,in position (1,2),we would look for neighbours in the positions (1,2),(1,1) and (0,1) in the down-oriented triangles array. Therefore, the natural way to store the boundary data is to use the diagonal and also to add a column before the first one and an extra row after the last one that we can see by applying the same idea to node 7, see Figure 2.28. Figure 2.28: A triangle structure stored in two arrays, one for up-oriented triangles and another for down-oriented triangles, with boundary data storage, B. Finally, we have obtained a data storage which does not differentiate from boundary nodes to internal ones and keeps the logical structure of a triangle. It has the problem that one needs to store the data by pairs (left side of the diagonal and right side of the diagonal), and that one will always work with two different arrays. Nevertheless, these two drawbacks are not a big problem in the solution of PDEs. The information is usually required by pairs: solution and right hand side, residual and previous step time solution. Regarding the use of two arrays, it may affect the performance if we are not cautious. Thus, we will always try to work first with one kind of triangles and later with the other one. 2.5.2 Extra-relaxation process in semi-structured grids The idea of applying multigrid on semi-structured grids is very simple. We will treat differently all the triangles that form the initial unstructured triangulation, see Figure 2.29. Nonetheless, we have to communicate them. The natural way is to consider each triangle to have Dirichlet boundary conditions, where these boundary conditions contain the overlap from neighbouring triangles, see Figure 2.30. However, 45 2. SCALAR MULTIGRID METHODS ON SEMI-STRUCTURED VORONOI GRIDS these “boundaries” have to be updated after every modification of their values, see Figure 2.31, for example, after the smoothing process and the prolongation update. Nevertheless, this methodology some times yields some problems that will require the use of an extra-relaxation process. (a) (b) Figure 2.29: (a) Initial unstructured grid of a domain. (b) Semi-structured grid obtained after one refinement level. Figure 2.30: Three isolated triangles after a virtual splitting. The first reason of considering an extra-relaxation process can be easily seen when considering the red-black smoother. To perform this method, we subdivide all the triangles in the grid, and therefore, the nodes into two different types: upand downoriented triangles, each of which can be updated simultaneously since each one only uses information from the other one. This means that, if we have two triangles of the initial grid sharing a common edge, and we want to apply a red-black smoother 46 2.5 Implementation Figure 2.31: Overlapping update process. we should keep the coherence across the triangles, see Figure 2.32 (a). Nevertheless, this can be very difficult if we consider not two triangles but hundreds of them. The splitting of the nodes for the red-black smoother applied on two coarse triangles is depicted in Figure 2.32 (b). We can see that the nodes at the overlapping regions belong to the same group (Black nodes). This means that one is not carrying out a true red-black smoother in the whole domain since the rule that defines a red-black smoother is not fulfilled by all the nodes. Hence, we are performing an “almost” redblack smoother, which, in the best case scenario will be as good as the true red-black smoother. Nonetheless, it seems that for the considered cell-centred discretizations we obtain a worse performance. (a) (b) Figure 2.32: (a) Red-Black smoother in two triangles keeping the coherence between them. (b) Practical red-black smoother. The second reason comes from the error smoothness that we obtain in the whole domain after smoothing each triangle by its own. We obtain a smooth error on the 47 2. SCALAR MULTIGRID METHODS ON SEMI-STRUCTURED VORONOI GRIDS triangles, but it seems not to be smooth enough at the nodes close to the common edges between the triangles. This effect is due to the fact that, as each triangle is treated separately, we obtain a smooth error on each single domain, but there can be a jump in the values from one triangle to its neighbours, which is sufficient to affect to the multigrid convergence factor. This behaviour can be easily appreciated in Figure 2.33, where we have depicted the error after two smoothing steps. The error on the whole domain seems smooth, but at the common edges of the triangles non smooth lines appear, deteriorating the convergence rate from a theoretical value of 0.10 to about 0.90.However, this effect can be even more dramatic, making the multigrid even to diverge in more complicated geometries. As conclusion, this issue must be solved in order to obtain an optimal method. Figure 2.33: Error after two smoothing steps, following the relaxation strategy presented in Figure 2.15. In order to get global smooth error it seems a good idea to apply a strong smoother in the nodes of the triangles that are close to the internal boundaries. In our case, the strongest considered relaxation scheme is the wormy smoother. Nevertheless, we would like to simultaneously update all the nodes, from both triangles, that are close to their common edge. Thus, we can consider a variant of the wormy smoother which smooths the lines at both sides at the same time, see Figure 2.34. However, we obtain a block-tridiagonal system of equations to solve, which makes this method really costly. Nonetheless, since the amount of nodes to be smoothed by this extra-relaxation process is of the order of √Totalnodes, the required amount of time to do that is small compared 48 3.1 Multigrid methods with discontinuous coefficients rates due to the weakness of the coarse-grid correction. Anyway, from the practical point of view, we will see that, in combination with strong smoothers, we will obtain very efficient multigrid methods. 3.1.1 Discretization on semi-structured grids of a diffusion problem with discontinuous coefficients We deal with the solution of the following diffusion equation: −∇·(κ(x, y)∇v) = f, in Ω,(3.1) v=g, on Γ,(3.2) where Ω is an open bounded domain, Γ is its boundary, fis a source term, ggives the Dirichlet boundary condition and κ(x, y) is the diffusion coefficient, that may be discontinuous. In particular, here we are interested in problems where κis piecewise constant. As done in Chapter 2, here we consider here a cell-centred finite volume discretization of problem (3.1)-(3.2) on a semi-structured grid. Firstly, we construct an acute Delaunay triangulation Ton the domain Ω,as the coarsest unstructured grid. We assume this triangulation to be fairly coarse, resolving the large-scale features of the domain, that is, this triangulation is taken fitting the geometry of the domain and the discontinuities of the diffusion coefficient in the way that on each triangle the value of κ remains constant. Secondly, we apply a regular refinement process to the triangles of T, that is, every triangle is divided into four congruent ones by connecting the midpoints of their edges, and so on until getting the desired fine scale to approximate the solution of the problem. Now, we are going to obtain the equation on the Voronoi point xccorresponding to a triangle Tof the triangulation T,see Figure 2.2. Note that Tis connected with three different triangles at most, T1, T2and T3.We denote by κcthe diffusion coefficient on the present triangle T, and by κithose coefficients on Ti.Using the divergence theorem on equation (3.1) at triangle T, we obtain −Z∂T (κ(x, y)∇v)·ndl =ZT f(x)dx, where nis the unit outward normal vector to ∂T and ∂T =l1∪l2∪l3.Then, by splitting the integral on the left-hand side into the sum of line integrals corresponding to edges 55 3. EXTENSION TO MORE COMPLEX MODELS l1, l2and l3,and by approximating them as the length of the corresponding edge, li, multiplied by the flux evaluated in the midpoint of the edge, we obtain the left-hand side of the equation. Regarding the right-hand side integral, it is replaced by the area of the triangle multiplied by the value of fat the Voronoi point. Finally, the equation reads: −1 meas(T) 3 X i=1 κH imeas(li)vh(xi)−vh(xc) di=f(xc),(3.3) where x1,x2,x3are the Voronoi points of the triangles adjacent to Twith common edges l1, l2and l3,respectively, and dithe distance between points xcand xi,with i= 1,2,3,(see Figure 2.2). The coefficients κH iappearing in (3.3) are the harmonic average between the corresponding diffusion coefficients, given by: κH i=2κcκi κc+κi ,(3.4) which is the most accurate method of known techniques of averaging (46, 56). Next, we are interested in obtaining the equations for the Voronoi points on the structured patches of the semi-structured grid. These equations follow from (3.3) by properly choosing the coefficients and using a double indexing for the notation of Voronoi points, more appropriate when structured grids are dealt with. Following the notation introduced in Section 2.1.2, given an arbitrary pair (k1, k2) associated with a node of Gh,the equations corresponding to the two Voronoi points xu k1,k2and xv k1,k2are given by Luu hvu h(xu k1,k2) + Lud hvd h(xd k1,k2) = fu h(xu k1,k2),(3.5) Ldu hvu h(xu k1,k2) + Ldd hvd h(xd k1,k2) = fd h(xd k1,k2),(3.6) where operators Luu h, Lud h, Ldu hand Ldd hcan be obtained from equation (3.3), and are given in stencil form as: 56 3.1 Multigrid methods with discontinuous coefficients Luu h=1 meas(T)       0 0 0 0 3 X i=1 κH i li di 0 0 0 0       , Lud h=1 meas(T)        0 0 0 0−κH 1 l1 d1−κH 3 l3 d3 0−κH 2 l2 d2 0        , Ldu h=1 meas(T)       0−κH 2 l2 d2 0 −κH 3 l3 d3−κH 1 l1 d1 0 0 0 0       , Ldd h=1 meas(T)       0 0 0 0 3 X i=1 κH i li di 0 0 0 0       , where the distances d1, d2, d3and the lengths l1, l2, l3are defined depending on the orientation of the triangle, as seen in Figure 3.1. For example, for an up-oriented triangle d2is defined as the distance between xu k1,k2and xd k1,k2−1, and l2as the length of the edge between those Voronoi points. (a) (b) Figure 3.1: Notation used to construct the stencil on a Voronoi point at (a) an up-oriented triangle or at (b) a down-oriented triangle. 3.1.2 Block-wise Multigrid For this type of problems, the same methodology considered in the previous chapter will be used. That is, a block-wise multigrid in which different smoothers are chosen for the triangles of the coarsest grid, is considered. Due to the fact that a different diffusion coefficient does not affect the anisotropies of the grid, we will follow the guideline depicted in 2.15 for the choice of the smoothers. Regarding the coarse-grid correction, as in the previous chapter, simple and local inter-grid transfer operators 57 3. EXTENSION TO MORE COMPLEX MODELS have been chosen in order to reduce to the bare minimum the communication between the coarsest triangles. In particular, the natural injection is selected as the prolongation operator Ik k−1; in fact, this is the only possible operator that is local for cell-centred discretizations, and the restriction is chosen as the adjoint of the prolongation. In this way, only the four fine-grid points surrounding a coarse-grid point are related to it by both inter-grid transfer operators. We would like to remark that the choice of these operators preserves four-point stencils when Galerkin approach is applied, which is very convenient for the simplicity of the connections between the triangular blocks of the grid. When large jumps in the diffusion coefficient κoccur in the domain, a direct discretization on coarse grids may not work properly (49, 56). However, in most of the multigrid methods proposed for discontinuous coefficient problems, the Galerkin approach has provided satisfactory results. This means that the coarse-grid operator is defined in terms of the fine-grid operator, Lk,the restriction, Ik−1 kand the prolongation, Ik k−1,in the following way: Lk−1=Ik−1 kLkIk k−1.(3.7) However, when simple injection is selected for the inter-grid transfer, the Galerkin operator (3.7) results to be inconsistent with the differential operator, and it must be replaced by the following expression: Lk−1=1 2Ik−1 kLkIk k−1.(3.8) This modified Galerkin approach provides better multigrid convergence rates than the original Galerkin operator (3.7), see (32, 40, 59). 3.1.3 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 the modified Galerkin approach. In the first experiment, problem (3.1) is solved on the unit square with two different distributions of diffusion coefficients. In the second one, the same problem is solved on a more complex domain. In the first numerical experiment, the proposed multigrid algorithm is applied by using an F-cycle with two preand two post-smoothing steps, 58 3.1 Multigrid methods with discontinuous coefficients whereas in the second one the influence of the type of cycle and the number of preand post-smoothing steps on the performance of the multigrid algorithm will be investigated. The computer used for the timing results in the subsequent experiments is an AMD at x2 at 2.9 GHz. 3.1.3.1 Diffusion problem on the unit square with discontinuous coefficients In the first experiment, we solve problem (3.1)-(3.2) on the unit square with different diffusion coefficients defined piecewise on two subdomains of different shapes, as shown in Figures 3.2(a) and 3.3(a). It consists of a benchmark problem taken from the literature (56). More concretely, in the first test case, the inner subdomain has a rhombus shape with a 0.5 side length, see Figure 3.2(a); whereas the second one consists of an hexadecagon with diameter 0.5 simulating a circle, see Figure 3.3(a). The right-hand side is defined as f(x, y) = xy, the Dirichlet boundary conditions are given by g(x, y) = x2+y2,and the diffusion coefficient values are κ= 0.333 ×105for the internal subdomains and κ= 2 for the rest of the domain, see Figures 3.2(a) and 3.3(a). In the same figures the corresponding coarsest grids are also represented. (a) (b) Figure 3.2: (a) Coarsest unstructured mesh for the first test case, and distribution of diffusion coefficients: κ= 0.333 ×105at the yellow region and κ= 2 at the green part. (b) Different smoothers for the triangles of the coarsest grid: white corresponds to red-black smoother, diamond smoother is represented by red, and wormy smoother by blue. 59 3. EXTENSION TO MORE COMPLEX MODELS (a) (b) Figure 3.3: (a) Coarsest unstructured mesh for the second test case, and distribution of diffusion coefficients: κ= 0.333 ×105at the yellow region and κ= 2 at the green part. (b) Different smoothers for the triangles of the coarsest grid: white corresponds to red-black smoother, diamond smoother is represented by red, and wormy smoother by blue. 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 3.2(b) and 3.3(b). Regarding the obtained multigrid convergence, in Table 3.1, the number of iterations necessary to reduce the initial residual in a factor of 10−10 are shown for both test cases. We observe an hindependent convergence for both problems, and although these results are slightly worse than those obtained in the case of constant diffusion coefficients. As expected, 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. 3.1.3.2 Diffusion problem on a composite material In the second experiment, problem (3.1)-(3.2) 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 3.4(a). The considered coarsest grid is shown in the same figure, and also we can observe that it is composed of triangles with very disparate shapes. For this reason, different smoothers are considered for the different triangles of the coarsest triangulation. In particular the smoothers chosen for these 60 3.1 Multigrid methods with discontinuous coefficients Levels Unknowns Rhombus Circle 4 6912 8 8 5 27648 8 9 6 110592 8 9 7 442368 8 9 8 1769472 9 9 Table 3.1: Number of iterations to reduce the initial residual in a factor of 10−10 for both test cases: the inner rhombus subdomain (Figure 3.2) and the inner circle subdomain (Figure 3.3). triangles are shown in Figure 3.4(b). In this way, the proposed block-wise multigrid is used for solving this problem. (a) (b) Figure 3.4: (a) Coarsest unstructured mesh and distribution of diffusion coefficients for the second experiment. Yellow color represents κ= 0.001 and green κ= 1. (b) Different smoothers for the triangles of the coarsest grid: white corresponds to red-black, diamond smoother is represented by red and wormy smoother by blue. Firstly, we want to compare the behaviour of the multigrid algorithm by considering both, direct discretization on coarse grids and the Galerkin approach. For this purpose, in Figure 3.5, the history of the convergence of the method by using an F(2,2)-cycle, for different numbers of refinement levels, is displayed. The stopping criterion is chosen as the final maximum residual to be less than 10−7.We can see that 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 reaches the desired value. Note, that in this experiment the use of Galerkin coarse-grid operator becomes mandatory. Next, we want to investigate the performance of the block-wise multigrid method 61 3. EXTENSION TO MORE COMPLEX MODELS Figure 3.5: Comparison between direct discretization and Galerkin approach for different numbers of refinement levels. proposed here, depending on the type of cycle and the number of smoothing steps. Varying these parameters, the number of iterations necessary to reduce the initial residual in a factor of 10−10 and the CPU time, together with the asymptotic convergence factor are shown in Table 3.2, for different number of refinement levels. In the case of an F-cycle, it is observed that the convergence is independent of the space discretization parameter, and that in few iterations the residual is reduced as desired. Regarding the behaviour of the V-cycle, it is clearly observed a very quickly deterioration in the asymptotic convergence factor with respect to the size of the mesh, even achieving divergence in some cases. This degradation of the convergence in the case of a V-cycle may be due to the lack of accuracy of the chosen inter-grid transfer operators. Notice that the price of obtaining an h-independent convergence rate is the use of an F-cycle. On the other hand, very simple inter-grid operators are used facilitating the communication between two connected patches, as well as strong and cheap smoothers are carefully chosen depending on the shape of the triangle. 3.2 Dominant convection The combination of diffusion and convection processes is extensively used. For instance, we can find many examples related to flow in porous media, including oil recovery, tracking of contaminants in groundwater flows and nuclear waste storage. Many schemes are known to fail if the convection is dominant, yielding solutions which are physically meaningless. The main difficulty comes from the interplay of the convection, a large scale effect, and diffusion, that comes from the Brownian motion which is a small scale 62 3.2 Dominant convection 5 levels 6 levels 7 levels 8 levels it (ρ) CPU it (ρ) CPU it (ρ) CPU it (ρ) CPU 2V 31(0.71) 17.29 131(0.95) 167.39 - - - - F 16(0.30) 19.59 16(0.30) 39.51 16(0.32) 96.99 17(0.32) 306.64 3V 13(0.58) 10.03 21(0.78) 36.90 - - - - F 11(0.17) 17.98 11(0.18) 36.77 11(0.19) 90.11 11(0.18) 277.02 4V 9(0.48) 8.88 16(0.65) 35.70 35(0.83) 217.04 - - F 8(0.15) 16.36 9(0.14) 37.99 9(0.15) 92.85 9(0.15) 273.17 Table 3.2: Number of iterations to reduce the initial residual in a factor of 10−10,corresponding asymptotic convergence rates (between brackets) and CPU-time for different numbers of refinement levels, by using a V-cycle and an F-cycle, and different numbers of smoothing steps ν= 2,3,4.The symbol −indicates divergence of the method. effect. Among all the numerical methods to stabilize the convection-diffusion problems we will use the simplest one, which is the upwind discretization. An overview of multigrid methods for dominant convection problems can be seen in (19, 49), where the authors point out some inherent problems that multigrid has when the convection is dominant, and also how to overcome those difficulties. 3.2.1 Discretization on unstructured triangular grids In the numerical experiment presented in Section 2.4.2, we presented a convectiondiffusion model, discretized by a central difference scheme, which was solved for a dominant diffusion situation. Nevertheless, it is more common to face a physical phenomena where the convection is dominant over the diffusion. In this case, the discretization used in Section 2.4.2 requires the use of very fine grids to obtain a physically acceptable solution, making its resolution prohibitive (49). On the other hand, the central difference scheme is a second order discretization, which is a desirable property for a scheme. However, in practice, this scheme is not always possible to use. When the convection is dominant, the discretization must be coherent with the direction of the flow. Using central differences, we consider all the directions to have the same importance, which is not true when there is a velocity field. As an alternative, we use an upwind discretization for the convective term, that despite having only first order, it will make the method stable no matter how big is the convective term. The equation for the convection-diffusion model written in divergence form is: 63 3. EXTENSION TO MORE COMPLEX MODELS ∇·(∇v−bv) = f, in Ω,(3.9) where b(x) is a given velocity field, whose divergence is assumed to be zero. We will use the same finite volume strategy used in the previous chapter. By using the divergence theorem on a triangle Tof the triangulation T,the following balance equation holds: Zl1 (∇v−bv)·n1dl1+Zl2 (∇v−bv)·n2dl2+Zl3 (∇v−bv)·n3dl3=ZT f(x) dx.(3.10) where niis the unit outward normal vector to the corresponding edge liof triangle T. Each of these line integrals can be split in two parts, the diffusion one, and the convective term. The diffusion part can be approximated as we did in Section 2.1.1. Regarding the convective part, we have the following line integral: −Zli bv·nidli,(3.11) which can be approximated by the length of the edge multiplied by an approximation of bv·niat the midpoint of the edge i. In this case, to approximate the value of v, an average between the values stored in the two triangles that share that common edge can be considered. This approach yields a second order approximation. However, it will not work when convection is dominant. In order to perform an upwind discretization, we have to approximate vby using only the value from one of those two triangles, following the direction of the flow. We can define a general formula for the upwind discretization: −bi(ξ vh(xc) + (1 −ξ)vh(xi)) meas(li),ξ= 1,if bi>0 ξ= 0,if bi≤0,(3.12) where biis the evaluation of b·niat the midpoint of edge i,xcis the node of the triangle in which we are currently discretizing, and xiis the node of the neighbouring triangle touching the corresponding edge, and ξis a parameter which is one or zero depending on the direction of the velocity. If we assemble the whole equation, we obtain: 1 meas(T) 3 X i=1 meas(li)vh(xi)−vh(xc) di−bi(ξ vh(xc) + (1 −ξ)vh(xi))=f(xc), with ξdefined as: (ξ= 1,if bi>0 ξ= 0,if bi≤0, (3.13) 64 3.3 Time dependent convection-diffusion problem b=0and = 1.We use a random initial guess and a zero right hand side to avoid round-off errors. Regarding the rest of the multigrid components, we have used the same inter-grid transfer operators presented in Section 2.2.1, direct discretization on coarse grids and two preand two post-smoothing steps. Since we want to compare the efficiency of the method, we will perform the experiments with only one time step. The time discretization parameter is fixed as τ= 0.01.The spatial grid is composed of 16384 triangles. In Table 3.6, we can see the results for an almost-right triangular grid and for a triangular grid characterized by a small angle. In all the cases the results are better than those obtained in Section 2.2.3. Despite the fact that the improvement is small, it is enough to realise that the convergence rate now also depends on the time discretization step. More concretely, it depends on the relation between τand the space discretization parameter. From this and considering the fact that the time discretization step can only improve the performance of the multigrid method, we can conclude that the strategy presented in Figure 2.15 to choose suitable smoothers is also valid for the time dependent convection-diffusion equation, since the steady problem would be the worst case. Right-Triangle Sharp triangle Smoother F(2,2) V(2,2) F(2,2) V(2,2) Red-Black 0.55 – 0.93 0.94 w-Red-Black 0.43 – 0.91 0.91 Diamond 3.6×10−23.7×10−20.87 0.88 Wormy 6.3×10−26.6×10−22.7×10−22.7×10−2 Table 3.6: Asymptotic convergence rates for an almost right triangular grid characterized by angles (45, 85), and for an isosceles triangular grid characterized by angles (85,85), by using four different smoothers. The time discretization is τ= 0.01 and a grid with 16384 nodes is considered. Now, we are going to consider an equilateral triangular grid. The same experiment as presented in Table 3.3 is performed for the time dependent convection-diffusion equation, with τ= 0.01.The corresponding results are shown in Table 3.7, and one can see that in all the cases the results have improved. 71 3. EXTENSION TO MORE COMPLEX MODELS b= (0,1) b= (0,−1) b= (1,1) b= (−1,−1) Green diamond 4.3×10−20.10 0.13 9.5×10−2 Red diamond 0.19 5.3×10−80.12 2,7×10−2 Black diamond 4.2×10−20.10 4.6×10−80.22 Green wormy 4.2×10−21.3×10−81.2×10−82.9×10−2 Red wormy 1.7×10−80.10 1.2×10−89.3 10−2 Black wormy 4.2×10−21.3×10−80.13 2.0×10−8 Red-Black 0.10 0.20 0.24 0.13 Table 3.7: Convergence factors for different smoothers using a V(1,1)-cycle with different flow directions and τ= 0.01. 3.3.3 Multigrid results on semi-structured grids Finally we are going to show the suitability of the proposed multigrid when semistructured triangular grids are considered. To this purpose, we solve equation (3.16) in a polygonal domain with three isolated holes inside. We have considered the circulatory flow given in equation (3.15). For these experiments we will fix = 1.The coarsest mesh is formed by 108 triangles with different geometries, and it is shown in Figure 3.7. The choice of the smoothers has been done by using the guide shown in Figure 2.15, and they are displayed in the same figure. In order to test the effects of the time dependent term, we will consider two different time discretization steps. In Table 3.8 we show the results for τ= 0.01 and in Table 3.9 for τ= 0.00001.We present the number of iterations to reduce the initial residual in ten orders of magnitude, the necessary CPUtime by using an AMD X2 at 2.9 GHz and the asymptotic convergence rates. When τ= 0.01,F-cycle presents an h-independent convergence rate, while that of V-cycle deteriorates as we increase the number of nodes. Despite this effect, the application of V-cycle seems to be more interesting as it requires less computational cost. Regarding the case in which τ= 0.00001,we observe that the convergence rate for Vand F-cycles is pretty similar. Thus, the use of the F-cycle seems completely unnecessary for small time steps. 72 3.3 Time dependent convection-diffusion problem (a) (b) Figure 3.7: (a) Grid after two refinement levels (b) Different smoothers for the triangles of the coarsest grid: white corresponds to red-black, diamond smoother is represented by red and wormy smoother by blue. V(2,2) F(2,2) Levels Unknowns Cycles ρhTime (s) Cycles ρhTime (s) 4 27648 10 0.12 1.7 9 0.10 1.7 5 110592 10 0.13 3.8 9 0.10 4.5 6 442368 11 0.15 11.8 9 0.11 13.7 7 1769472 11 0.17 39.5 9 0.12 44.7 8 7077888 11 0.18 190.8 9 0.12 221.4 Table 3.8: Number of iterations to reduce the initial residual in a factor of 10−10,corresponding asymptotic convergence rates and CPU-times for different numbers of refinement levels, by using V-cycle and F-cycle, where τ= 0.01. 73 3. EXTENSION TO MORE COMPLEX MODELS V(2,2) F(2,2) Levels Unknowns Cycles ρhTime (s) Cycles ρhTime (s) 4 27648 4 6.3×10−40.6 4 6.3×10−40.8 5 110592 6 1.3×10−22.3 5 1.3×10−23.1 6 442368 8 5.3×10−28.4 8 5.3×10−212.1 7 1769472 9 0.11 32.2 9 0.10 45.2 8 7077888 10 0.12 161.2 9 0.11 207.6 Table 3.9: Number of iterations to reduce the initial residual in a factor of 10−10,corresponding asymptotic convergence rates and CPU-times for different numbers of refinement levels, by using V-cycle and F-cycle, where τ= 0.00001. 74 Chapter 4 Multigrid methods on staggered grids Until this chapter we have only considered scalar equations. Nevertheless, this is very limited, since in many areas the use of scalar fields is insufficient. Consider, for example, computational fluid dynamics where the velocity is represented by vectors. Moreover, we may be more interested in the interaction between the velocity of a fluid and other variables like the pressure, or the concentration of a contaminant. Therefore, in this chapter we are going to deal with the use of vectors and scalar values together. Besides, we will also have to deal with the anisotropy created by the meshes. For that, we will adapt the multigrid components that we have developed for scalar equations to work with staggered grids. 4.1 Discretization on unstructured triangular grids In this section, a discretization for the following system of equations is presented: u+OP=f O·u= 0 , in Ω,(4.1) u·n= 0, on ∂Ω,(4.2) where uis the velocity field, and Pa scalar value denoting the pressure. We will approximate the normal component of the velocity at the midpoint of the edges, see Figure 4.1, while the scalar variables will be approximated at the Voronoi points. For 75 4. MULTIGRID METHODS ON STAGGERED GRIDS the normal components of the velocity, we will use the same notation as in Section 3.2. Regarding the properties of the triangulation, we assume the same restrictions established in Section 2.1.1. Figure 4.1: Unstructured mesh and its associated Voronoi grid with dots for the scalar values, and squares denoting the places where the projection of the vector values are stored. For the continuity equation, we will use the discretization technique called VAGO, presented in (50). In this way, the resulting discretization for any triangle Treads as: 1 meas(T) 3 X i=1 (meas(li)ui)=0,(4.3) where uiis the numerical approximation of the outward normal component of the velocity on edge i. Once we have discretized the scalar equation, now we deal with Darcy’s law, whose discretization at the midpoint of edge iof triangle Tis: ui+P(xi)−P(xc) di =fi.(4.4) Here, uiis the normal projection of the velocity on edge i,P(xc) is the pressure at the Voronoi point of triangle T, P(xi) the pressure corresponding to the Voronoi point xi, and fiis the normal projection of the right hand side evaluated in the middle point of the edge i. As seen in (50),the resulting scheme is a first order discretization. 76 4.2 Multigrid method 4.2 Multigrid method 4.2.1 Coarse-grid correction In this subsection we are going to deal with the prolongation and restriction for discretizations based on edges. For cell-centred nodes we have used the same transfer operators introduced in Section 2.2.1. In order to define the restriction and the prolongation for discretizations based on edges, we are going to distinguish three different type of nodes, one for each edge. Then the restriction and prolongation for each type of node only involves grid-points of its same type. 4.2.1.1 Restriction We considered a six-point restriction operator with weights of 1 4or 1 12 depending on the relative location of the nodes, see Figure 4.2 (a). For example, for the grid-points located at edges of type L1: rc(xic,jc) = (3rf(xi,j)+3rf(xi,j+1)+3rf(xi+1,j+1) + rf(xi,j−1) +rf(xi+1,j) + rf(xi+2,j+1))/12,(4.5) where rcand rfare the residual grid functions defined on the coarse and fine grid respectively, and the coarse indices icand jcare related to those on the fine grid by the following expressions: i= 2ic−1 and j= 2jc−1. 4.2.1.2 Prolongation Regarding the prolongation, given a fine-grid node belonging to an up-oriented coarse triangle, we apply a simple injection from the coarse-grid point located at the same coarse triangle. In other case, we apply an average of the three closest coarse-grid points, see Figure 4.2 (b). 4.2.2 Smoothers Since the considered system is a saddle point problem, it is well known that standard smoothers do not provide good convergence, and then more sophisticated smoothers have to be considered. More concretely, Vanka type smoothers (51), as well as distributive smoothers (58), are common choices for these type of problems. Next, we develop novel Vanka type smoothers adapted to staggered discretizations on structured 77 4. MULTIGRID METHODS ON STAGGERED GRIDS (a) (b) Figure 4.2: (a) Restriction for a L1 type node (circles), where the red nodes have a weight of 1 4and the green ones of 1 12 . (b) Prolongation for a L1 type node, where some nodes are updated by injection, whereas others, like the orange node, need to be calculated from the three coarse grid surrounding nodes. triangular grids. Due to the properties of the discrete operators considered here, we can also design a suitable distributive smoother. 4.2.2.1 Vanka type smoothers In general, Vanka type smoothers consist of simultaneously updating a set of unknowns. Of course, these blocks can be updated in different orderings. As a generalization of the smoothers introduced in Chapter 2, next we are going to present three different Vanka type smoothers. Vanka red-black smoother: For this smoother a block of four unknowns belonging to the same triangle is updated. In particular, the three velocity unknowns located at the edges, denoted here as v1, v2, v3,and the scalar unknown associated with the Voronoi point, represented by P, are simultaneously updated. 78 4.2 Multigrid method Following the idea of a red-black ordering, we first update all the blocks corresponding to up-oriented triangles, and secondly all the blocks associated with down-oriented triangles. Notice, that the unknowns located at the edges are relaxed twice per smoothing step, see Figure 4.3. Figure 4.3: Vanka red-black smoother for staggered triangular grids. The application of this smoother requires us to solve for each triangle the following small 4 ×4 system of equations:        1001 dv1 0101 dv2 0011 dv3 l1 meas(T) l2 meas(T) l3 meas(T)0            δv1 δv2 δv3 δP    =    r1 r2 r3 rs    ,(4.6) where riis the residual on the node corresponding to the edge iand rsthe residual corresponding to the scalar unknown. Vanka diamond smoother: In a Vanka diamond smoother we have to update simultaneously all the unknowns located at two neighbouring triangles. In this case, we have to solve a 7 ×7 system of equations for each block: 79 4. MULTIGRID METHODS ON STAGGERED GRIDS                   1001 dv1 0 0 0 0101 dv2 0 0 −1 dv2 0011 dv3 0 0 0 l1 meas(T) l2 meas(T) l3 meas(T)0 0 0 0 0 0 0 0 1 0 1 dv5 0 0 0 0 0 1 1 dv4 0−l2 meas(T)0 0 l5 meas(T) l4 meas(T)0                               δv1 δv2 δv3 δP1 δv4 δv5 δP2             =             r1 r2 r3 rs1 r4 r5 rs2             . (4.7) Vanka wormy smoother: In this case, we have to update simultaneously the unknowns located at a row of triangles. This means that, whereas a tridiagonal system had to be solved in the scalar case, now a block tridiagonal system must be solved, which means that the computational cost will be very high. 4.2.2.2 Distributive smoother In this kind of smoother, the idea is to decouple the system of equations by introducing a right preconditioner in the smoothing procedure, in order to apply fast smoothers to the decoupled system, see Figure 4.4. More concretely it consists of two steps: •Predictor step: After decoupling the system, we smooth the new variables by smoothing the scalar and vector unknowns separately. •Corrector step: The final step consists in updating the original variables by using the values obtained in the predictor step. Next, we are going to describe the concrete distributive smoother for our system: IhGh Dh0! Uh Ph!= fh 0!.(4.8) By introducing the new variables Uh Ph!= IhGh 0−Ih! U0 h P0 h!,(4.9) the new decoupled system yields: 80 4.4 A numerical experiment on semi-structured grids Figure 4.11: Different smoothers for the triangles of the coarsest grid: diamond smoother is represented by red, and wormy smoother by blue. In Table 4.3 results similar to those presented in the previous experiment are shown. Again very good results are obtained in this case, which makes this strategy also suitable in practice. Finally, as we expected from the study performed on the previous section, the strategy based only on distributive relaxation is prefered due to its lower computational cost. Therefore, this methodology will be used in next chapter for solving the most real problem considered in this thesis. 87 4. MULTIGRID METHODS ON STAGGERED GRIDS V(2,2) F(2,2) Levels Unknowns Cycles ρhTime (s) Cycles ρhTime (s) 3 24192 10 0.12 18.3 9 8.9×10−225.6 4 96768 10 0.13 45.2 9 8.9×10−261.8 5 387072 11 0.12 148.6 10 9.0×10−2174.3 6 1548288 11 0.12 491.7 10 9.8×10−2528.4 7 6193152 11 0.12 1710.4 11 9.7×10−21879.1 Table 4.2: Number of cycles, asymptotic convergence factors and CPU times for different type of cycles and different numbers of refinement levels, by using the global distributive smoother. V(2,2) F(2,2) Levels Unknowns Cycles ρhTime (s) Cycles ρhTime (s) 3 24192 14 0.25 25.7 11 0.16 31.6 4 96768 14 0.23 80.6 12 0.16 96.8 5 387072 14 0.23 207.5 12 0.16 203.5 6 1548288 14 0.23 718.3 12 0.14 727.7 7 6193152 14 0.22 2399.7 13 0.16 2512.2 Table 4.3: Number of cycles, asymptotic convergence factors and CPU times for different type of cycles and different numbers of refinement levels, by using the smoothers depicted in Figure 4.11. 88 Chapter 5 Multigrid method on semi-structured grids for the resolution of the Darcy-Oberbeck-Boussinesq model 5.1 Density driven flows Density driven flows are flows carried out by differences in the density of the fluid (27). These flows may appear in all the situations in which there is a variation of density or when there are two immiscible fluids with a different density. Also, if there is only one phase, density variations across the fluid may be produced by a gradient of temperatures or a concentration of a solute. These types of flows appear frequently in porous media, where the speed of the stream is characterized by Darcy’s law. More concretely, when underwater bears are used to store hot water in summer to be pumped out in winter. The introduced water, that is in steady state, has a different temperature and therefore, a gradient of temperatures will appear creating a slow flow driven by the different densities. Another example of density driven flow in porous media appears in the storage of a chemical waste underground. The presence of the chemical material in the underwater bear creates a flow due to the variation of its concentration, which 89 5. MULTIGRID METHOD ON SEMI-STRUCTURED GRIDS FOR THE RESOLUTION OF THE DARCY-OBERBECK-BOUSSINESQ MODEL may lead to a leakage of the residual. Clean coal technologies are another example of density driven flows in porous media. This technology consists in introducing the CO2 produced by the combustion in underwater saline aquifers. Nevertheless, not in all the cases where there is a variation of the density a density driven flow models the behaviour of the process. To consider a density driven flow, the distribution of the density must be coupled with the velocity of the flow. As the density is directly related with the temperature or with the concentration of a solute, then we consider this variable the variable of interest. 5.1.1 The Darcy-Oberbeck-Boussinesq model In order to model density driven flows, we are going to work out the Darcy-OberbeckBoussinesq model, DOB from now on, by using the basic equations to model a fluid. We start with Navier-Stokes equations: ρDv Dt =−∇P+µ∇2v+ρgz,(5.1) ∂ρ ∂t +∇·(ρv)=0,(5.2) where ρis the density of the fluid, Pis the pressure, gthe gravity force, µthe viscosity, vthe velocity field, zis the vector (0,−1) in Cartesian coordinates, and D Dt denotes the material derivative. Now, we consider the usual assumptions to achieve the Stokes equation, incompressibility, creeping flow and also stationary. Then, the material derivative of the left side of equation (5.1) and the derivative of ρare equal to zero, and then we obtain the Stokes equations: −µ∇2v+∇P=ρgz,(5.3) ∇·v= 0.(5.4) We can consider the fluid to have an almost constant volume, which means that assuming the divergence of the velocity to be zero, condition (5.4), is acceptable even though we are going to consider variations of density later. In (8), it is shown that the variation of the volume is very small for the majority of the fluids, below 10−3. Regarding the viscosity term, if we consider an isotropic porous media, an approximation of this expression can be done by considering a linear response of the viscous term: µ∇2v=−µφv κ,(5.5) 90 5.1 Density driven flows in which φis the porosity, and κthe permeability. In order to obtain the “Darcy velocity”, we multiply the porosity by the velocity field: u=φv.(5.6) Introducing equations (5.5) and (5.6) into equation (5.3), we obtain the Darcy’s law: u=−κ µ(∇P−ρgz).(5.7) Darcy’s law was formulated by Henry Darcy after some experiments about the flow of water through sand, and it was published in (10). The change of the density of a fluid can be produced by a variation of the temperature, pressure or concentration of a solution. We will consider only the variation of concentration as our intention is to model the sequestration of CO2.Therefore, the next equation to consider is the advection equation: DC Dt =D∇2C. (5.8) In this equation, Dis the diffusion coefficient relative to the concentration C. Next, we proceed to expand the material derivative: ∂C ∂t +v·∇C=D∇2C. (5.9) Substituting vin equation (5.9), by using equation (5.6), we obtain the desired equation: φ∂C ∂t =−u·∇C+φD∇2C. (5.10) Despite the fact that we have considered an incompressible fluid, we must take into account the variations of density in the buoyancy term carried out by the variation of concentration of a solute. Then, the variation of density follows the equation: ρ=ρ0+δρC, (5.11) where the density ρis a linear function of concentration in which ρ0is the density of the lighter fluid and δρ is the difference between the density of the two fluids. 91 5. MULTIGRID METHOD ON SEMI-STRUCTURED GRIDS FOR THE RESOLUTION OF THE DARCY-OBERBECK-BOUSSINESQ MODEL 5.2 The Darcy-Oberbeck-Boussinesq model in CO2storage We are going to focus on the use of the DOB model in the simulation of storage of CO2 in underwater saline aquifers. Hence, the variation of the density of the flow will be due to the dissolution of a solute, CO2, in the water. Firstly, we want to explain the clean coal technology, so we can understand its interest. Next, we will adapt the DOB model to this particular case. 5.2.1 The CO2sequestration process The carbon capture and storage is a process that consists of obtaining the CO2from a source, usually from an important source like a power plant. After that the CO2is compressed, transported and finally stored. 5.2.1.1 Carbon capture There are three capturing processes: Post-combustion, pre-combustion and oxy-fuel. The first, post-combustion, is the simplest one. This method burns the fuel in the old-fashion way. After the combustion, there is a process to capture the CO2from the flue gases. This method is very interesting because it is an external procedure to the burning process, which means that the old power stations can be upgraded with this technology to turn them into clean power stations. Pre-combustion consists in capturing the CO2before the fuel is burned. In this process, initially the fuel is partially oxidized, turning the fuel into carbon monoxide, hydrogen and water. The separation before the combustion yields cleaner gases to deal with, and also the gases can be compressed, making this last process simpler. However, this method cannot be applied to old power stations. The third option, oxy-fuel combustion, consists in burning the fuel directly with oxygen instead of air. The resulting flue gases only contain CO2and water. Therefore, all the gases can be treated after the process, not only a fraction of the total flue gases, as it was the case in the first method. 5.2.1.2 Compression and transport Once the CO2is obtained from the flue gases, it must be compressed to a super-critical state to facilitate its transportation. The compression means that part of the energy 92 5.2 The Darcy-Oberbeck-Boussinesq model in CO2storage obtained in the power station will have to be spent in that process. The amount of necessary energy is about the 40% of the energy produced. Hence, the cost of the electricity will be increased in an important factor. However, it is believed that this percentage can be reduced to 20%. After the compression of CO2, it has to be transported. For this part of the process, many existing pipelines, which are already in use for other purposes, can also be used to pump CO2. 5.2.1.3 Storage mechanisms The final step of the sequestration process is to store the CO2.For that there are many possibilities. One possibility is to dissolve CO2in the oceans, however, the introduction of CO2in the water changes its acidity. We could consider storing it in tanks, however, this would lead to a situation close to the storage of nuclear waste, since these residuals, despite not being so dangerous are long-living as well. Also, there is always the possibility of these tanks to explode. A better option which has been used since many years, is to store CO2in the reservoirs where the oil is stored. Introducing CO2in the reservoirs started in order to facilitate the extraction of oil. Despite the fact that this method could be considered as an option, the main problem is that not all countries have oil and hence such reservoirs to store the CO2. Another option is to store CO2underground, in other kind of places. Among all the possible places, the most interesting is the use of saline aquifers. Saline aquifers have been already used to store industrial wastes, since they are of no use for humankind and also they are very extended, making easy the transportation. The use of saline aquifers is also very interesting because the use of water to dilute CO2means that after the necessary time the CO2will turn into a mineral, making the storage permanent. However, before this final storage, the CO2have to be saved for several years. In that time, the CO2must be stored without leakages to the surface. To understand better the effectiveness of the storage method, we are going to see in more detail the way in which the CO2is stored. Firstly, the CO2is pumped into the saline aquifer. Here, there are four mechanisms that avoid the CO2to return to the surface: •Structural trapping: The CO2may try to scape to the surface, however, it is stored due to the walls of the formation. 93 5. MULTIGRID METHOD ON SEMI-STRUCTURED GRIDS FOR THE RESOLUTION OF THE DARCY-OBERBECK-BOUSSINESQ MODEL •Dissolution trapping: The CO2can be dissolved in the water. When this happens, the density of the water is changed making the water rich in CO2to go to the bottom. •Residual trapping: This trapping occurs when pumping the CO2. The porous rock acts like a sponge, then the supercritical CO2displaces the fluid stored in the rock. However, after that initial displacement, the fluid returns, displacing the majority of the CO2, but leaving behind some residual CO2stored in the pores of the rock. •Mineral trapping: This trapping process is the slowest one, but it is permanent. When the CO2is dissolved into water, it creates carbonic acid. After a long time, this carbonic acid reacts with the rocks of the boundaries to form a mineral, finally trapping the CO2. 5.2.2 Nondimensionalization of the Darcy-Oberbeck-Boussinesq model We present the nondimensionalization of the Darcy-Oberbeck-Boussinesq model: u=−κ µ(∇P−ρgz),(5.12) φ∂C ∂t =−u·OC+φD∆C, (5.13) ∇·u= 0,(5.14) ρ=ρ0+δρC. (5.15) For the nondimensionalization of this model, we are going to follow the process presented in (24). Firstly, we need to define the basic parameters to be used in the nondimensionalization : µ, u,H (the domain width), κ, φ and δρ. After selecting the basic parameters, we define the dimensionless parameters: •u∗=κ[m2]δρ[Kg m3]g[m s2] µ[Kg s m ] •P∗=µ[Kg s m ]u[m s]H[m] κ[m2]=δρgH •t∗=φ[−]H[m] u[m s] 94 5.2 The Darcy-Oberbeck-Boussinesq model in CO2storage •ρ∗=δρ[Kg m3] Using these parameters we can define new variables: Pn=P P∗,un=u u∗, tn=t t∗, xn=x H, zn=z H.(5.16) The final step is to introduce the nondiomensional variables into the original system. We begin with equation (5.12): unu∗=−κ µ(∇n H(PnδρgH)−ρgz),(5.17) and after some algebraic manipulations we obtain the following equation: un=−κ µu∗δρg(∇n H(PnH)−ρ0+δρC δρ z).(5.18) As δρg =µu∗ κwe can modify the equation as follows: un=−(∇nPn−ρ0z δρ −Cz).(5.19) Finally, if we regroup all the terms affected by the gradient in a new variable P0= Pn−ρ0zn δρH ,where we have introduced the variable znin the gradient in order to keep the equation the same, we obtain the following dimensionless equation: un=−(∇nP0−Cz).(5.20) Next step is the nondimensionalization of the equation (5.13), by substituting the old variables by the dimensionless ones: φ∂C ∂tn u∗ φH =−unu∗∇nC H+φD∆n H2C. (5.21) Multiplying the equation by H/u∗we can remove the coefficients of almost all terms: ∂C ∂tn=−un∇nC+φD u∗H∆nC. (5.22) Equation (5.22) can be simplified to contain a single parameter, since φD u∗His equal to the inverse of the Rayleigh parameter. Therefore, the final equation reads: ∂C ∂tn=−un∇nC+1 Ra∆nC. (5.23) 95 5. MULTIGRID METHOD ON SEMI-STRUCTURED GRIDS FOR THE RESOLUTION OF THE DARCY-OBERBECK-BOUSSINESQ MODEL It is very easy to get the dimensionless equivalent of equation (5.14): ∇n H·un u∗=0=∇n·un.(5.24) If we finally assemble all the equations, we obtain the desired dimensionless equation system: ∂C ∂tn=−un∇nC+1 Ra∆nC, (5.25) un=−(∇nP0−Cz),(5.26) ∇n·un= 0.(5.27) 5.3 Solver strategy In this section, we explain the numerical strategy that we have used to solve the DOB model. Regarding the discretization methods used to discretize the equations of the model, we have used the techniques explained in previous chapters. Regarding the non-linearity in the term −un∇nC, we use a fixed point method. Thus, we will split the system into two sub-systems: ∂C ∂tn=−un old∇nC+1 Ra∆nC, (5.28) and un=−(∇nP0−Coldz) (5.29) ∇n·un= 0.(5.30) The fixed point method follows the next scheme, also represented in Figure 5.1: •Using the multigrid method explained in Section 3.3, the parabolic equation (5.28) is solved using an old guess of un.In this way, a new approximation of Cis calculated. •Next, using the already calculated Cwe solve the system of equations (5.29) and (5.30) with the multigrid method that we presented in Section 4. •After obtaining a new velocity, we will repeat the process until the concentration does not change its value more than 1% from one step to the next one. 96 5.4 Numerical experiments (a) (b) Figure 5.9: CO2concentration after 9.97 dimensionless seconds for the two considered meshes. Rayleigh number equal to 10000 In this experiment, we increase the Rayleigh number in one order of magnitude. Hence, we expect a result driven more by the convection than by the diffusion process. The mesh used is shown in Figure 5.2 (a). In this case, as we reduce the diffusion term, we must be careful with the artificial diffusion introduced by the discretization method. If it becomes bigger or similar to the real diffusion, then the obtained results will not be acceptable. In order to know whether the artificial diffusion may modify our results or not we study the order of magnitude of the artificial and real diffusion. As it is known, the introduced artificial diffusion is proportional to the space discretization parameter. In our case, the initial unstructured triangulation contains 24 triangles, and each one contains 4096 scalar nodes. This means that the total amount of scalar nodes is about 105.As the Rayleigh number is 104and the domain measures are 0.5 high and 1 width, then we can say that the space discretization is about 10−5.Therefore, we can conclude that the artificial diffusion is of one order of magnitude smaller than the real diffusion. Now, we proceed to study the behaviour of the experiment. In this case, after 0.07 dimensionless seconds, we can see in Figure 5.10 (a) that some fingers appear in the boundaries and also in the centre of the domain. Therefore, we can conclude that the instabilities have started earlier than for the previous case. Besides, more fingers are created, see Figure 5.10 (b). 103 5. MULTIGRID METHOD ON SEMI-STRUCTURED GRIDS FOR THE RESOLUTION OF THE DARCY-OBERBECK-BOUSSINESQ MODEL (a) (b) Figure 5.10: (a) CO2concentration after 0.07 dimensionless seconds. (b) CO2concentration after 0.25 dimensionless seconds. In this experiment, the bigger fingers do not ingest the small ones so easily. Thus, after 0.39 dimensionless seconds, we can see in Figure 5.11 (a) that the majority of the fingers are growing. Nonetheless, some have been absorbed reducing the initial number of fingers, see Figure 5.11 (b) and Figure 5.12 (a). In the last picture, Figure 5.12 (b), some fingers have already reached the bottom. Also we can observe that some fingers are almost detached from the upper boundary. Nevertheless, the connection is still kept, although it is weak. (a) (b) Figure 5.11: (a) CO2concentration after 0.39 dimensionless seconds. (b) CO2concentration after 0.68 dimensionless seconds. 104 5.4 Numerical experiments (a) (b) Figure 5.12: (a) CO2concentration after 1.22 dimensionless seconds. (b) CO2concentration after 2.50 dimensionless seconds. 5.4.1.2 Pipe shape domain Finally, we perform two experiments in a pipe-shape domain in order to test the influence of the shape of the domain in the creation of fingers. In the first experiment, we impose a low Rayleigh value, making the behaviour of the test to be driven by the diffusivity. In the second experiment, a Rayleigh value of 5000 is considered. A domain with a curve part on the boundary is considered in order to see how the shape of the domain influences the CO2absorption process. In Figure 5.13 the smoothers used for the resolution of the advection equation are depicted. Figure 5.13: Different smoothers for the triangles of the coarsest grid: white corresponds to red-black, diamond smoother is represented by red and wormy smoother by blue. Rayleigh number equal to 100 In this example, we expect a high diffusion component since the Rayleigh number is very 105 5. MULTIGRID METHOD ON SEMI-STRUCTURED GRIDS FOR THE RESOLUTION OF THE DARCY-OBERBECK-BOUSSINESQ MODEL low and the instabilities appear for a Rayleigh number higher than 380 (27). Looking at Figures 5.14, 5.15 and 5.16, we can see that in this case, a single finger has been created. (a) (b) Figure 5.14: (a) CO2concentration after 0.03 dimensionless seconds. (b) CO2concentration after 1.00 dimensionless seconds. (a) (b) Figure 5.15: (a) CO2concentration after 2.08 dimensionless seconds. (b) CO2concentration after 2.51 dimensionless seconds. 106 5.4 Numerical experiments (a) (b) Figure 5.16: (a) CO2concentration after 5.02 dimensionless seconds. (b) CO2concentration after 9.99 dimensionless seconds. Rayleigh number equal to 5000 Again, we face a high Rayleigh value. As expected, several fingers appear, see Figure 5.17. Later, they start to grow and merge between themselves, see Figure 5.18. Finally, due to the shape of the domain a large finger is created that ingest all the new fingers that have been created, see Figure 5.19. (a) (b) Figure 5.17: (a) CO2concentration after 1.02 dimensionless seconds. (b) CO2concentration after 2.08 dimensionless seconds. 107 5. MULTIGRID METHOD ON SEMI-STRUCTURED GRIDS FOR THE RESOLUTION OF THE DARCY-OBERBECK-BOUSSINESQ MODEL (a) (b) Figure 5.18: (a) CO2concentration after 3.00 dimensionless seconds. (b) CO2concentration after 4.02 dimensionless seconds. (a) (b) Figure 5.19: (a) CO2concentration after 6.07 dimensionless seconds. (b) CO2concentration after 8.02 dimensionless seconds. 5.4.2 Numerical experiments in heterogeneous media This section deals with domains where different layers of soil are considered. We perform two experiments in which we will study the effect of these layers in the absorption of CO2. 5.4.2.1 Hourglass shape domain In this first experiment, we consider an hourglass shape domain in which the layer at the middle of the hourglass has a Rayleigh number value of 8000, while the rest of the domain has a value of 1000, see Figure 5.20. The considered smoothers for the resolution of the advection equation are depicted in Figure 5.20 (b). 108 5.4 Numerical experiments (a) (b) Figure 5.20: (a) Rayleigh number values in the considered domain. Blue stands for a Rayleigh number of 1000 and red for a number of 8000. (b) Different smoothers for the scalar case: white corresponds to red-black, diamond smoother is represented by red and wormy smoother by blue. Before the fingers reach the layer at the centre of the domain, the behaviour of the experiment is similar to the homogeneous media experiments, see Figure 5.21. Nevertheless, once the fingers are approaching the centre layer, they start to reduce their speed,see Figure 5.22 (a), and even they retract themselves until a bigger finger is created, see Figure 5.22 (b). Finally, once a finger of CO2has reached the other section of the hourglass, Figure 5.23 (a), the stream of CO2is stable and a single finger of CO2 goes through the rest of the domain, Figure 5.23 (b). (a) (b) Figure 5.21: (a) CO2concentration after 1.09 dimensionless seconds. (b) CO2concentration after 2.47 dimensionless seconds. 109 5. MULTIGRID METHOD ON SEMI-STRUCTURED GRIDS FOR THE RESOLUTION OF THE DARCY-OBERBECK-BOUSSINESQ MODEL (a) (b) Figure 5.22: (a) CO2concentration after 3.56 dimensionless seconds. (b) CO2concentration after 5.01 dimensionless seconds. (a) (b) Figure 5.23: (a) CO2concentration after 7.02 dimensionless seconds. (b) CO2concentration after 9.99 dimensionless seconds. 5.4.2.2 Strata domain Now, we want to analyse the behaviour of the process when many different layers are considered in the domain. The jump in the Rayleigh number between the layers is different in all the cases. Figure 5.24 (a) shows the four different layers. The light green area has a Rayleigh value of 100 (diffusive response), the blue one has a value of 1000, the red zone a value 5000 and the yellow of 8000. For this experiment, we have used 1024 scalar nodes per coarse triangle, and a time step of 0.0001,while the rest of the parameters are the same as before. For the advection equation, the smoothers used are depicted in Figure 5.24 (b). In Figure 5.25 (a), we see an initial stable layer of CO2. However, once the CO2reaches the zone of a higher Rayleigh value, Figure 5.25 (b), 110 5.4 Numerical experiments (a) (b) Figure 5.24: (a) Rayleigh numbers in the considered domain. Light green represents a Rayleigh number of 100, Blue stands for 1000, red for 5000 and yellow for 8000. (b) Different smoothers for the scalar case: white corresponds to red-black, diamond smoother is represented by red and wormy smoother by blue. many fingers start to grow, which continue growing, as we seen in Figure 5.26 (a). In this case, when the fingers reach the layer with a Rayleigh number value of 5000, the fingers continue their movements without any visible effect, Figure 5.26 (b). In order to study whether the layer with a Rayleigh number value of 5000 actually modifies the results or not, we decided to run this experiment with only two layers, that with a Rayleigh value of 100 was kept unchanged while the other three layers were merged into a single one with a Rayleigh number value of 1000. Comparing those two experiments, some differences in the behaviour of the fingers appeared. In the case of four layers, the fingers tended to reduce their number, but on the other hand they were bigger. In Figure 5.27 (a) and (b) the fingers have reached the bottom of the domain; some plumes of low concentration start to appear close to the boundaries. Finally, in Figure 5.28 (a) and (b) we can see that as the domain is almost full, the dissolution of CO2becomes slower. 111 5. MULTIGRID METHOD ON SEMI-STRUCTURED GRIDS FOR THE RESOLUTION OF THE DARCY-OBERBECK-BOUSSINESQ MODEL (a) (b) Figure 5.25: (a) CO2concentration after 0.33 dimensionless seconds. (b) CO2concentration after 0.83 dimensionless seconds. (a) (b) Figure 5.26: (a) CO2concentration after 1.17 dimensionless seconds. (b) CO2concentration after 2.33 dimensionless seconds. (a) (b) Figure 5.27: (a) CO2concentration after 3.33 dimensionless seconds. (b) CO2concentration after 4.99 dimensionless seconds. 112 En relaci´on al material presentado en los Cap´ıtulos 4 y 5, nuestra intenci´on es publicar un art´ıculo tras la presentaci´on de esta tesis. 119 6. CONCLUSIONS 120 Appendices 121 Appendix A Basic multigrid The efficient numerical solution of many real problems modelled by partial differential equations strongly depends on the resolution of the corresponding large linear systems resulting from their discretization. Since their development in the 60’s, multigrid methods have been proved to be among the most efficient numerical algorithms for solving this kind of systems, achieving asymptotically optimal complexity at least for elliptic problems. That is, the required computational work to solve a discrete problem is of the order of the number of unknowns of the corresponding system. These methods are iterative solvers with a characteristic property: the convergence of a multigrid method is independent of the discretization grid size. Multigrid methods are mainly based on the acceleration of the convergence of common iterative methods by using solutions obtained on coarser meshes as corrections. Most basic relaxation schemes suffer when smooth components are present in the error, however these components become high frequency components on a coarser grid, so that it makes sense to go down to such a grid and relax directly on the error by using the residual equation. We then return to the fine grid to correct the approximation obtained there with the correction from the coarser grid. Two principles are involved in the development of multigrid methods: the first one is the fact that some classical iterative methods have a strong smoothing effect on the components of the error corresponding to the high frequencies (high oscillating error components). This effect is due to the fact that relaxation schemes (as iterative methods are known in the multigrid context) basically average neighbouring discrete values, leading to a reduction of the oscillating character of the error. However, when the error is mainly composed of low frequency components, the differences between 123 A. BASIC MULTIGRID neighbouring discrete points are small and this leads to a slowdown in the convergence of these methods (7). In Figure A.1 this behaviour is shown, by choosing a random initial guess and applying some iterations of a classical iterative method, in particular, a Gauss-Seidel method. It is observed that in a few iterations the oscillating character of the algebraic error, is reduced resulting in a very smooth error. Notice that although the error is smooth, it is not significantly reduced in magnitude. The components asInitial error Error after 10 iterations Figure A.1: Gauss-Seidel smoothing effect on the error. sociated with the low frequencies are annihilated by relaxation methods in a very slow way. When the error becomes smooth, the convergence of these methods gets worse. In order to avoid this, the second fundamental idea of these methods is employed, which is that a smooth error can be well represented on a coarser grid, where all the computations can be done faster because of the reduction of the number of grid-points. This strategy can be applied in a recursive way using coarser and coarser grids until a level of refinement is reached on which the application of an exact solver would be feasible. Consequently, multigrid methods are based on using a hierarchy of grids in order to eliminate all the different components of the error. There are two possibilities to achieve an efficient interplay between the smoothing and the coarse–grid correction, which give rise to two different approaches to multigrid algorithms: geometric and algebraic. Algebraic multigrid (AMG) appears suitable for problems on unstructured grids. The use of semi-structured grids appears as an alternative to the use of AMG for unstructured grids for relatively complex domains, since these grids allow the efficient implementation of a geometric multigrid algorithm. 124 A.1 Multigrid components A.1 Multigrid components Defining a sequence of grids G0, G1, . . . , Gl,characterized by grid sizes h0> h1> . . . > hl,the multigrid iteration operator, Mk,can be recursively defined as follows: M0= 0, Mk= (Sk)ν2Ck(Sk)ν1, k = 1, . . . , l, where (A.1) Ck=Ik−Ik k−1(Ik−1−(Mk−1)γ)L−1 k−1Ik−1 kLk, where Skis the smoother on grid Gk, Ckrepresents the coarse-grid correction part of the algorithm, and the fact that the coarse grid equation is approximately solved by applying γiterations of a multigrid cycle, is indicated by the term (Ik−1−(Mk−1)γ)L−1 k−1(49). For defining (A.1), a discrete operator Lk:G(Gk)→G(Gk) approximating Ll,and a suitable relaxation process Sk:G(Gk)→G(Gk) must be defined on each grid k= 0, . . . , l −1,and the corresponding transfer grid operators between consecutive grids are denoted by Ik+1 k:G(Gk)→G(Gk+1),(A.2) Ik k+1 :G(Gk+1)→G(Gk).(A.3) The idea of using a sequence of grids permits us to process them in different ways. These manners of processing the grids are determined by the cycle index γ, which indicates the number of multigrid steps to perform on coarser grids. The value γ= 1 yields the V-cycle, which is the easiest recursive definition of a multigrid cycle. By increasing the number of times coarser grids are processed, other type of cycles can be defined. For example, γ= 2 yields the so-called W-cycle, and another type of cycle, as the F-cycle, is also commonly used. The F-cycle consists of going down to the coarsest grid and recursively interpolating to the next finer grid and applying a V-cycle until the finest mesh is reached (49). The different types of cycles are represented in Figure A.2. From the definition of the multigrid operator (A.1), it is deduced that many details are open for discussion and decision. All its components, that is, the smoother, intergrid transfer operators (restriction and prolongation), the choice of the coarser grids and the operators on them, the type of cycle and the number of preand post-smoothing 125 A. BASIC MULTIGRID Figure A.2: Types of cycles. steps, have to be properly chosen. For each concrete problem they have to be specified, and although for some problems the components which give good results are known, when dealing with new applications it is always a challenge to find the most suitable ones. The coarsening strategy defines the hierarchy of grids necessary for the implementation of the geometric multigrid algorithm. Different ways of building the coarser grids can be considered, the more standard being the strategy of doubling the mesh size from Gkto Gk−1,that is, hk−1= 2hk. Once the hierarchy of grids has been built, suitable discrete operators Lkon each coarse grid Gk, k = 0, . . . , l −1 have to be chosen approximating the fine-grid discrete operator Ll.The most natural way to define these operators is to use the direct discretization of the equations on each coarse grid, and it is known as discretization coarse-grid approximation (DCA). This choice usually gives rise to reasonable approximations to Ll,ensuring the overall consistency of the discrete problem. An alternative choice of Lkis the so-called Galerkin coarse-grid approximation (GCA), which consists of using the inter-grid transfer operators, Ik+1 kand Ik k+1,to define Lkfrom Lk+1,that is, Lk=Ik k+1Lk+1Ik+1 k.(A.4) When Galerkin coarse-grid operators are used, the restriction and the prolongation are often selected to be adjoint to each other for symmetric positive definite operators. This approximation is necessary for problems with discontinuous coefficients, and also in the application of algebraic multigrid, for example. The choice of inter-grid transfer operators is, of course, closely related to the coarsening strategy. Transfer of information between two consecutive grids Gk,and Gk+1 is done by the restriction and prolongation operators, Ik k+1 and Ik+1 k,which map 126 A.2 Local Fourier Analysis (k+1)−grid functions on k−grid functions and k−grid functions on (k+1)−grid functions, respectively. The simplest way to define a prolongation operator is through linear interpolation. Other possible choices for the prolongation can be chosen, like high-order interpolations, for example. More sophisticated interpolation methods are sometimes helpful, depending on the problem and on the quantity to be transferred. For instance, matrixdependent prolongations (13, 14), whose construction is based on the discrete operator under consideration, are useful for problems with discontinuous coefficients and in the algebraic multigrid framework. On the other hand, the simplest restriction operator which comes into mind is the injection operator. Injection is easily programmed, but restriction operators which define a coarse-grid value by a certain averaging of neighbouring fine-grid values provide robust algorithms for more general problems. The full-weighting restriction is commonly used. The smoother usually plays an important role in multigrid algorithms, above all in the geometric approach. From the early stages of multigrid development until now, a wide variety of smoothing procedures has been investigated apart from the standard ones like Jacobi and Gauss-Seidel type relaxations in all their varieties. Smoothers as ADIor CG-type methods, the incomplete LU-decomposition (ILU), distributive relaxations, Vanka-type smoothers, and so on, have been widely used in the multigrid framework, and this list could be extended with many more, (49). A.2 Local Fourier Analysis The local Fourier analysis (LFA), introduced by Brandt in 1977 (4), is a quantitative analysis for multigrid algorithms, in the sense that it provides accurate predictions of the asymptotic convergence rates of such methods. This analysis is based on the Fourier transform theory, and a good introduction can be found in the books by Trottenberg et al. (49), and Wienands and Joppich (57). A k-level local Fourier analysis is a very useful tool to predict the convergence rate of a multigrid algorithm. The main idea of this analysis is formally to extend all multigrid components to an infinite grid, neglecting the boundary conditions, and to restrict the analysis to discrete linear operators with constant coefficients. Despite these restrictions, if boundary con127 A. BASIC MULTIGRID ditions are appropriately treated, in general this analysis matches the numerical results satisfactorily. In this Appendix, our intention is just to explain the easiest application of LFA, that is, the smoothing analysis. This analysis is enough to study the performance of the smoothers and to understand the LFA concept. Nevertheless, it is not enough to simulate the multigrid behaviour. For that, at least a two-grid analysis has to be carried out. For a deeper understanding of this analysis we refer the reader to (57), where LFA is explained in detail. A.2.1 Smoothing analysis In order to explain the smoothing analysis, we will consider the standard finite difference discretization of Laplace operator on a Cartesian grid. LFA is based on the idea that a discrete grid function can be written as a formal linear combination of Fourier modes. ϕ(θ,x) = eιθ·x/h,(A.5) where xis a point of the grid and θis a frequency in (−π, π]2. It is well-known that for each frequency the corresponding Fourier mode is an eigenfunction of the discrete Laplace operator, with corresponding eigenvalue, called symbol in this context, ˜ Lh(θ) = 1 h2(4 −eιθ1/h −eιθ2/h −e−ιθ1/h −e−ιθ2/h).(A.6) Any standard smoothing operator can be described by a splitting of the discrete operator Lh=L+ h+L− h.Then, the symbol of any relaxation operators is given by ˜ Sh(θ) = −˜ L− h(θ) ˜ L+ h(θ).(A.7) Since the idea of smoothing analysis is to study how the smoother eliminates the high-frequency components of the error, we will only consider ˜ Sh(θ), for θbeing a high frequency, see Figure A.3. Then, the smoothing factor is defined as the maximum of the modulus of these values. For example, Gauss-Seidel method provides a smoothing factor of 0.5 which is obtained for the frequency (θ1, θ2)=(π/2,0.6435). The smoothing analysis is a good start to know a priori the behaviour of a method. 128