SciELO - Scientific Electronic Library Online

 
vol.45 número1Photocatalytic activity of ZnAl2O4 spinel for procion red degradation under UV irradiationA hybrid group based re-key management scheme for secure communication in wireless sensor índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Latin American applied research

versión On-line ISSN 1851-8796

Lat. Am. appl. res. vol.45 no.1 Bahía Blanca ene. 2015

 

A multigrid method for the solution of composite mesh problems

S.S. Sarraf, E.J. Lopez, G.A. Rios Rodriguez and V.E. Sonzogni

Departamento de Mecánica Aplicada, Facultad de Ingeniería, Universidad Nacional del Comahue, CONICET, Buenos Aires 1400, Q8300IBX Neuquén, Argentina. sssarraf@gmail.com, ezequiel.lopez@fain.uncoma.edu.ar
Centro de Investigación de Métodos Computacionales (CIMEC), Universidad Nacional del Litoral, CONICET, Predio Conicet-Santa Fe, Colectora Ruta Nac. 168 / Paraje El Pozo, S3000GLN Santa Fe, Argentina. gusadrr@yahoo.com.ar, sonzogni@intec.unl.edu.ar

Abstract— The Composite Finite Element Mesh method is useful for the estimation of the discretization error and, in addition, for the nodal solution improvement with a small increase in the computational cost. The technique uses two meshes with different element size to discretize a given problem and, then, it redefines the resulting linear system. On the other hand, Multigrid methods solve a linear system using systems of several sizes resulting from a hierarchy of meshes. This feature motivates the study of the application of the Multigrid strategy together with the Composite Mesh technique. In this work, it is proposed a Multigrid method to solve problems where the Composite Mesh is applied. The goal of the proposal is to achieve both, the advantages of the Multigrid algorithm efficiency and the solution improvement given by the Composite Mesh technique. The new method is tested with some elliptic problems with analytical solution.

Keywords— Multigrid Method; Composite Mesh; Numerical Solution Improvement.

I. INTRODUCTION

In this work, a Multigrid (MG) technique able to solve the linear system arising from the application of the Composite Mesh (CM) strategy is presented. The application of the Finite Element Method (FEM) to discretize a partial differential equation leads to a linear system, which could be solved with some MG strategy in a very efficient way. Of course, the element size used to discretize the problem domain directly affects both the computational cost, through the size of the linear system, and the error level. For a given problem, the idea behind the CM strategy is to perform a linear combination between the discrete systems of equations computed with FEM using two meshes of different element size. These meshes must have nodes in common. Typically, one mesh is the homogeneous refinement of the other one. With an appropriate choice of the coefficients of the linear combination, the CM technique could give a better nodal solution than those obtained from each mesh individually without increasing the computational cost significantly (Bergallo et al., 2000). Although the CM method has been tested for several kind of problems (Sarraf, 2011), the best performance of the strategy was found for elliptic problems where the solution has a high degree of regularity (Sonzogni et al., 1996; Bergallo et al., 2000).

The MG method is a well known strategy useful for solving differential equations using a hierarchy of meshes defined over the problem domain (Lallemand et al., 1992; Wesseling, 1992; Koobus et al., 1994; Venkatakrishnan and Mavriplis, 1994; Mavriplis, 1995; Chan et al., 1997; Briggs et al., 2000; Arnold, 2001; Okusanya, 2002; Kim et al., 2004). MG algorithms are applied to a wide range of problems, primarily to solve linear and nonlinear boundary value problems. One of the first applications of the MG techniques was to elliptic problems, which remain today as one of the typical applications of the method.

Since both MG and CM methods use meshes with different degree of refinement grouped in a hierarchy of grids and are suitable for elliptic problems, the motivation to integrate these strategies arises. In this paper we describe the proposed method, which is applied to some elliptic test problems with analytical solution on unstructured meshes, where discretization errors are analyzed.

II. MULTIGRID METHOD

Let Ω be a bounded domain in ℜd with boundary ∂Ω, d=1, 2, 3 being the space dimension. Consider the following elliptic differential equation with homogeneous boundary conditions

(1)

Assume that the operator L is self-adjoint, i.e.
(Lu, v)=(u, Lv) for any u, v H L2(Ω) and that it is positive in the sense that (Lu, u)>0 for all uH, u≠0, where the subspace H contains smooth functions which vanish on ∂Ω. With these properties, to solve the boundary value problem given by Eq. (1) is formally equivalent to minimize the quadratic functional

(2)

The problem can be rewritten in compact form as follows

(3)

which means to find the argument that minimizes F over all the functions in H.

Given a triangulation of Ω with element size h, denoted by Ωh, let Hh be the finite subspace of H consisting on the functions uh which are continuous in Ω, polynomial in each element and vanish on the boundary domain. Then, the discrete problem is written as follows

(4)

Since the operator L is positive and self-adjoint, the problem of determining the function uH which satisfies Eq. (3) and (Lu,v)=(f,v) for all vH are equivalent (Briggs et al., 2000). Hence, to solve Eq. (4) is equivalent to find uh Hh so that

for all (5)

In order to solve Eq. (5) in its weak form, let εhi be a function of a base for Hh, with εhi such that εhi(Nj)=δij, Nj being the nodes of Ωh. Then, choosing the test functions vh as the basis functions εhi and after assembling all the rows of the matrix and their corresponding right-hand side, a discrete system of the form Ahuh=fh is obtained.

Given two meshes, where the first is an homogeneous refinement of the second one, the two-grid iteration is the basis for building the MG method. In summary, this iteration is composed by the following steps:

  1. A few iterations with some iterative method such as Gauss-Seidel or damped Jacobi are performed in order to smooth out the residue . This step will reduce the high frequency components of the error but not completely the low frequency ones. The new residue, which seems relatively smooth, can be properly approximated on a coarser mesh.
  2. The residue is projected on the coarse mesh where the error equation is solved. This action minimizes the low frequency components of the error because they behave as high-frequency ones on the coarse mesh.
  3. Finally, this solution on the coarse mesh is transferred to the fine mesh where it is used to improve the approximation of step 1.

Firstly consider the relaxation step, where the main goal is to provide an inexpensive method for the elimination of oscillatory errors in the approximation . This objective can be achieved carrying out local changes in the following way: , where s∈ℜ is a suitable step size. The choice of s is carried out in the sense of minimizing the functional over all the possible choices, i.e.

(6)

The abstract formulation for the correction process on the coarse mesh can be stated as follows: Let HHHh be the space of the coarse grid, i.e. the set of piecewise polynomial functions associated with a standard coarse grid ΩH. The goal is to correct the approximation with a function that will approximate the new smooth error. The correction is : . The choice of is made in order to find the best correction in the sense of the functional minimization, i.e.

(7)

In the process of adding a function HH to a function Hh, appropriate coefficients must be found in order to write as a function of Hh. Then, the prolongation operator IHh is sought to carry over (Briggs et al., 2000).

To determine the coarse grid operator AH corresponding to the operator Ah, it is necessary to work with nodal vectors and to transfer the minimization principle to matrix terms. Let Ni, i=1,..., M be the nodes of the mesh Ωh. Then, the system Ahuh=fh is equivalent to the minimization problem

(8)

where . Similarly, the problem of the coarse mesh correction is equivalent to the matrix minimization principle

(9)

Therefore, the version of Ah in the coarse mesh is obtained by performing the computation AH=(IHh)T Ah IHh where (.)T denotes transpose matrix. The matrix (IHh)T takes a vector from the fine grid and gives a vector in the coarse mesh, i.e. it is a restriction operator. Consequently, it makes sense having IhH = (IHh)T.

Finally, let rh=fhAh and fH=IhHrh. Then it is obtained

(10)

Therefore, since is independent of wH, to minimize is equivalent to minimize FH(wH) on the vectors wH∈ℜM/2, which is the correction scheme for the coarse mesh.

The MG algorithm is obtained by applying the two grid iteration recursively. Specifically, the problem in the coarse mesh is solved by a two-grid iteration involving a coarser mesh. This procedure is applied at each mesh level until a coarse enough grid is obtained in which the system could be solved by a direct method. A MG iteration from the finest grid to the coarser one and returning to the fine one, is called cycle. The exact structure of a cycle depends on the value of γ, the number of iterations of two grids in each intermediate step. The case γ=1 is called V-cycle, while γ=2 is called W-cycle. The main procedure of the MG method is given by Algorithm 1, where in the conditional line 6, represents the Euclidean norm for vectors.

The computation of the mesh hierarchy and the matrices for each level is performed by Algorithm 2, where the container A is a vector in which the i-th "element" is the matrix of the i th MG level.

The pseudo-code for the MG solver is given by Algorithm 3. Here, function smoother applies some iterative method in order to smooth out the residue. The functions restrict and prolong perform the tasks of restricting and interpolating of vectors between grids of different levels, respectively.

III. COMPOSITE MESH TECHNIQUE

We consider elliptic problems with the following form

(11)

where Ω⊂ℜd and n is the outward unit normal to ∂Ω. In general, the boundary domain is formed by parts with Dirichlet conditions (ΓD) and Neumann conditions (ΓN=∂Ω\ΓD). The function f(x) is the source term and µ(x)> 0 ∀x∈Ω is the diffusivity coefficient.

The composite finite element mesh applied to elliptic problems can be used to improve the numerical solution without an appreciable increase in the computational cost and also to estimate the discretization error (Sonzogni et al., 1996; Bergallo et al., 2000; Sarraf et al., 2007). In the h version, the method consists in replace the discrete operator for a given mesh (the fine mesh), by a linear combination of the discrete operators corresponding to this mesh and a coarser mesh with nodes in common with the first one. In this case, the interpolation polynomials retain the same degree in both meshes. Then, assuming that the fine mesh is obtained from the homogeneous refinement of the other grid, the connection between both meshes is forced with the shared nodes. The participation factor of each mesh in the compound model, i.e. the coefficient in the linear combination between the meshes, is introduced in such a way to minimize the discretization error.

Let ΩH be a discretization of the problem domain Ω and Ωh the mesh obtained by a homogeneous refinement of ΩH. The element size of meshes ΩH and Ωh are H and h, respectively. After discretizing the problem (11) with FEM, the systems of equations Ahuh=fh and AHuH=fH are obtained with meshes Ωh and ΩH, respectively. Now, we define the discrete operator AHh as a representation of matrix AH into the space of the matrix Ah. A natural choice for the computation of that matrix is to use the prolongation operator presented in the last section, as AHh = IHh AH. This definition is quite general, since the election of the prolongation operator could lead to different methods with dissimilar performance. In this work, we use the injection operator (Briggs et al., 2000), which has proven to work fine for the CM method (Sonzogni et al., 1996; Bergallo et al., 2000). The modified matrix proposed in the CM technique ACM is computed as a linear combination of matrices Ah and AHh. Vectors fHh and fCM are defined in an analogous way.

Then, the approximate solution by the CM method uHh is obtained from the following system (Sonzogni et al., 1996)

(12)

When α=1, the solution in the fine mesh is recovered and with the solution tends to the corresponding solution in the coarse grid. The coefficient α depends on the regularity of the analytical solution of the problem (Bergallo et al., 2000; Sarraf, 2011).

The asymptotic error of the numerical approximation has the form (Sonzogni et al., 1996)

(13)

where h is the mesh size, C is a constant independent of h and q>p. Therefore, an extrapolation analysis of the error leads to the following estimation (Sonzogni et al., 1996; Bergallo et al., 2000; Toro et al., 2005)

(14)

The improvement introduced by the CM method with respect to the FEM solution is realized in the nodal values of the solution and, thus, it must be evaluated using a discrete norm of the error. Let Ni, i =1,... ,M be the nodes of the fine mesh and Hh the discrete space associated with Ωh. The interpolant πhu of u in the space Hh is defined as

(15)

Then, the solution of Eq. (12) is a better approximation to πhu in Hh than the FEM solution. This fact is verified in the tests presented in section V.

IV. MULTIGRID METHOD TO SOLVE COMPOSITE MESH PROBLEMS

As pointed out above, the main objective of this paper is the integration of the MG and CM techniques for solving elliptic problems. Suppose that there are n levels involved at each MG iteration and that the mesh corresponding to the j-th level is obtained from the (j−1)-th homogeneous refinement of a given initial mesh (corresponding to level 1). Now, we introduce m "mixing" levels, 1≤m<n, in which the standard linear operators of the MG method are replaced by the linear combination proposed in the CM strategy. In the k-th level of mix, 1≤km, the grids of levels nk+1 and nk take part in the mesh composition. In other words, for level k: AkCM = αAn−k+1+(1−α)(IhH)kAn−k replaces matrix An−k+1, where An−k+1 and An−k are the matrices corresponding to levels n−k+1 and nk, respectively; and (IhH)k is the prolongation operator from level nk to level nk+1. The right-hand side vector of the highest grid level is replaced by the linear combination proposed by the CM strategy, i.e. fn=α fn+(1−α)(IhH) fn-1. Figure 1 outlines the strategy.


Figure 1: Multigrid for Composite Mesh (V -cycle).

The proposed method uses the same MG main algorithm (Alg. 1) changing the procedure which computes the linear systems called at line 2. The pseudo-code for the computation of the modified matrices is presented in Algorithm 4.

With the same set of parameters (smoother, number of pre- and post-smoothing steps, etc.) the resolution of each level increases the computational cost if the grid mixture is introduced for such level. This occurs because the mixture of the meshes increases the bandwidth of the matrix for a given level. Therefore, the number of mixing levels m should be kept as low as possible while keeping the full error reduction attainable by the CM technique. As will be shown in the numerical examples, the finer the mesh, the higher the decrease in the nodal error produced by the CM strategy with respect to FEM. Hence, an increase in the levels of mixing is expected with that of the grid levels.

V. NUMERICAL TESTS

In order to show the effectiveness of the proposed MGCM (Multigrid for Composite Mesh) strategy, we will consider three different two-dimensional problems with analytical solution. For these tests, we use in the analysis the standard l norm and the Euclidean norm for vectors. The error e is computed as the difference uuh, where u and uh are the nodal values of the analytical and the numerical solutions, respectively.

A. Poisson problem with constant coefficients

Consider the following problem in a quadrangular domain

(16)

The source term is defined such that the analytical solution is given by the following expression

In this case, the domain is discretized by an unstructured grid with 82 triangular elements and 52 nodes. This grid is the coarsest mesh of the problem and the grid sequence is obtained by homogeneous refinement. The parameters selected for this problem consist in three steps of pre- and post-smoothing of damped Jacobi with relaxation parameter 0.7 (Briggs et al., 2000). The tolerance applied to the residue for the convergence of the linear system is 1 × 10−6. The participation factor for the mesh composition is obtained taking into account the regularity of the exact solution, giving α =4/3 (Bergallo et al., 2000; Sarraf, 2011).

The number of iterations to achieve convergence (it.), the number of grid levels (n) and a measure of the elapsed computational time are shown in Table 1 for the V -cycle and inTable 2 for the W -cycle. In the MGCM case, the number of mixing levels (m) is also included.

Table 1: Results of the V -cycle for the Poisson problem with constant coefficients.

Table 2: Results of the W -cycle for the Poisson problem with constant coefficients.

The increased computational time of the MGCM method is due to the fact that the matrix of the system is less sparse than in the "standard" MG case. Regarding the behavior of the proposed strategy in section IV for solving the CM problem with MG, we conclude that the MGCM method preserves the features of the original MG method. This assertion is based on the fact that the number of iterations needed to reach convergence, with and without "mixture", seems to remain constant. Figure 2 presents the Euclidean norm of the residue as a function of the number of iterations for the case with five grid levels for the V and W-cycles. As can be observed, the curves for the MG and MGCM methods seem to be superimposed with a slightly difference for the W-cycle.


Figure 2: Residual norm as a function of the number of iterations for the Poisson problem with constant coefficients and five grid levels.

For the V-cycle with six grid levels, two levels of mesh composition are used in order to achieve the full error reduction that the CM technique should give. Nevertheless, for the W-cycle only one mixing level was sufficient in all the solved cases. In Fig. 3 the nodal error as a function of the discretization size (h) for both norms is presented. In this figure, we refer with FEM and CM as the MG and MGCM solutions, respectively.


Figure 3: (a) Euclidean and (b) infinity norm of the nodal error as a function of the discretization size h for the Poisson problem with constant coefficients.

A noticeable feature is the error reduction rate of the CM technique.

B. Poisson problem with variable coefficients

Consider the following elliptic problem with variable coefficients in the unit square:

(17)

where µ(x,y) = 1+xy2 and u(x,y) = sin(5πx) cos(3πy). The mesh and the whole set of parameters are the same as in the previous test. Tables 3 and 4 show the results for the V-cycle and W-cycle, respectively. In Fig. 4 it can be observed the reduction of the error when the CM technique is used. In this case the convergence rates of the proposed MGCM strategy and the standard MG method are similar, in particular for one mixing level.

Table 3: Results of the V-cycle for the Poisson problem with variable coefficients.

Table 4: Results of the W-cycle for the Poisson problem with variable coefficients.


Figure 4: Residual norm as a function of the number of iterations for the Poisson problem with variable coefficients using five grid levels.

Figure 5 shows the Euclidean and infinity norms of the nodal error as a function of the discretization size. In this case, it is mandatory to use two mixing levels for the V-cycle with six grid levels in order to achieve the expected rate in the error reduction.


Figure 5: (a) Euclidean and (b) infinity norm of the nodal error as a function of the discretization size h for the Poisson problem with variable coefficients.

C. Laplace problem in a L-shaped domain

Let the Laplace problem in the L-shaped domain Ω = (-1,0) × (-1,1)∪(0,1)×(0,1) where the exact solution in polar coordinates is given by

(18)

The mesh of the first level is an unstructured grid with 82 elements and 52 nodes. Again, we use three steps of pre- and post-smoothing of damped Jacobi with factor 0.7. The participation factor for the mesh composition is obtained taking into account the regularity of the analytical solution, resulting α =25/3/(25/3-1) (Grisvard, 1986; Sarraf, 2011). The methods are assumed to reach the convergence when the Euclidean norm of the residue is less than 1 × 10−6.

The number of iterations to achieve convergence (it.), the number of grid levels (n) and the elapsed computational time are presented in Table 5 and 6 for the V-cycle and W-cycle, respectively. For the MGCM strategy, the number of mixing levels m is also indicated.

Table 5: Results of the V-cycle for the Laplace problem in a L-shaped domain.

Table 6: Results of the W-cycle for the Laplace problem in a L-shaped domain.

Figure 6 shows the norm of the nodal error as a function of the mesh size. Again, a reduction of the error norm can be reached when the MGCM technique is applied. In this case the relative reduction of the error between MG and MGCM is smaller than in the previous example because of the lower regularity in the analytical solution.


Figure 6: (a) Euclidean and (b) infinity norm of the nodal error as a function of the iscretization size h for the Laplace problem in a L-shaped domain.

VI. CONCLUSIONS

In this work we propose a Multigrid method able to solve the linear system arising from the application of the Composite Mesh technique. The strategy is based on the substitution of the linear operator defined by the MG method in some "mixing" levels. The substitute linear operator in a given mixing level is computed as a linear combination of the matrix in the current level and the matrix corresponding to the immediate inferior (coarser) level. The coefficient of the linear combination depends on the regularity of the analytical solution of the elliptic PDE (Partial Differential Equation) problem. The results of the solved problems show that the new method introduces a slight degradation of the Multigrid convergence, which is reflected in the elapsed time for the resolution of the linear system. The best performance of the proposed method was found for the W-cycle, particularly in cases where the analytical solution of the problem has high regularity. Nevertheless, a great reduction in the nodal error is obtained with the MGCM technique respect to the "standard" MG method which, in this study, corresponds to the numerical solution of the Finite Element Method.

ACKNOWLEDGEMENTS
This work has received financial support from Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET, Argentina), Universidad Nacional del Comahue (UNComa, Argentina), Universidad Nacional del Litoral (UNL, Argentina, grant CAI+D 2009-III-4-2), Agencia Nacional de Promoción Científica y Tecnológica (ANCyT, Argentina, grants PICT 2010-2492/16, PICT-PRH 2009-0147), and was performed with the Free Software Foundation/GNU-Project resources as GNU-Linux-OS, GNU-Octave, as well other open source resources such as Xfig and LATEX.

REFERENCES
1. Arnold, D.N., A Concise Introduction to Numerical Analysis, University of Minnesota, 195-206 (2001).
2. Bergallo, M.B., C.E. Neuman and V.E. Sonzogni, "Composite mesh concept based FEM error estimation and solution improvement," Comp. Meth. App. Mech. Engng., 188, 755-774 (2000).
3. Briggs, W.L., V.E. Hemson and S.F. McCormick, A Multigrid Tutorial, Society for Industrial and Applied Mathematics (2000).
4. Chan, T.F., S. Go and L. Zikatanov, "Lecture Notes on Multilevel Methods for Elliptic Problems on Unstructured Grids," VKI 28th Computational Fluid Dynamics, 1-76 (1997).
5. Grisvard, P., "Boundary value problems in plane polygons. Instructions for use," E.D.F. Bulletin de la Direction des Etudes et Recherches, Serie C Mathematiques, Informatique 1, 21-59 (1986).
6. Kim, J., K. Kang and J. Lowengrub, "Conservative multigrid methods for ternary Cahn-Hilliard systems," Comm. Math. Sci., 2, 53-77 (2004).
7. Koobus, B., M.H. Lallemand and A. Dervieux, "Unstructured volume agglomeration MG: solution of the Poisson equation," Int. J. Num. Meth. Fluids, 18, 27-42 (1994).
8. Lallemand, M.H., H. Steve and A. Dervieux, "Unstructured multigridding by volume agglomeration: current status," Comp. Fluids, 21, 397-433 (1992).
9. Mavriplis, D.J., "Multigrid techniques for unstructured meshes," Institute for Computer Applications in Science and Engineering Report, 95-27 (1995).
10. Okusanya, T., Algebraic Multigrid for Stabilized Finite Element Discretizations of the Navier Stokes Equations, PhD Thesis, Massachusetts Institute of Technology (2002).
11. Sarraf, S.S., M.B. Bergallo and V.E. Sonzogni, "Problemas elípticos resueltos mediante mallas compuestas aplicando métodos Multigrilla," Mecánica Computacional, XXVI, 696-710 (2007).
12. Sarraf, S.S., Análisis y generalización de la técnica de malla compuesta, PhD Thesis, Fac. de Ing. y Cs. Hídricas, Universidad Nacional del Litoral (2011).
13. Sonzogni, V.E., M.B. Bergallo and C.E. Neuman, "Uso de una malla compuesta para estimar errores de discretización y mejorar la solución en elementos finitos," Mecánica Computacional, XVI, 123-132 (1996).
14. Toro, S., V.E. Sonzogni and C.E. Neuman, "Elementos finitos de diferentes órdenes para problemas de elasticidad plana y mezcla de sus mallas," Mecánica Computacional, XXIV, 3171-3185 (2005).
15. Venkatakrishnan, V. and D. Mavriplis, "Agglomeration Multigrid for the 3D Euler equations," AIAA, paper 94-0069, 94-69 (1994).
16. Wesseling, P., An Introduction to Multigrid Methods, John Wiley & Sons (1992).

Received: November 26, 2012. Received by the Subject Editor: May 15, 2014.
Accepted: June 19, 2014. Recommended by Subject Editor: Walter Tuckart.