Performance of Geometric Multigrid Method for Two-Dimensional Burgers’ Equations with Non-Orthogonal, Structured Curvilinear Grids

: This paper seeks to develop an efficient multigrid algorithm for solving the Burgers problem with the use of non-orthogonal structured curvilinear grids in L-shaped geometry. For this, the differential equations were discretized by Finite Volume Method (FVM) with second-order approximation scheme and deferred correction. Moreover, the algebraic method and the differential method were used to generate the non-orthogonal structured curvilinear grids. Furthermore, the influence of some parameters of geometric multigrid method, as well as lexicographical Gauss–Seidel (Lex-GS), η -line Gauss–Seidel ( η -line-GS), Modified Strongly Implicit (MSI) and modified incomplete LU decomposition (MILU) solvers on the Central Processing Unit (CPU) time was investigated. Therefore, several parameters of multigrid method and solvers were tested for the problem, with the use of non-orthogonal structured curvilinear grids and multigrid method, resulting in an algorithm with the combination that achieved the best results and CPU time. The geometric multigrid method with Full Approximation Scheme (FAS), V-cycle and standard coarsening ratio for this problem were utilized. This article shows how to calculate the coordinates transformation metrics in the coarser grids. Results show that the MSI and MILU solvers are the most efficient. Moreover, the MSI solver is faster than MILU for both grids generators; and the solutions are more accurate for the Burgers problem with grids generated using elliptic equations.


Introduction
Many Engineering problems concern complex geometries, in which the use of Cartesian, cylindrical or spherical coordinate systems is not practical or appropriate, and discretization on structured curvilinear or unstructured grids is preferred.
Structured curvilinear grids are based on the mapping of the physical domain over a simple shape computational domain, for instance, a rectangle [1]. The use of such grids simplifies the application of boundary conditions, since the conditions existing on the boundaries of the physical plane are transferred exactly to the boundaries of the modified plane, and thus do not need approximation. With this option, a system of global equations is achieved and it is possible to write and discretize differential equations in this new system.
The discretization of mathematical models by the Finite Volume Method (FVM) [2] approximates a system of differential equations through a system of algebraic equations of the form where A represents the matrix of coefficients, φ is the vector of unknowns and b is the vector of independent terms. The FVM can be applied to any type of grid, in complex geometries and in different coordinate systems [3]. The attractive feature of this method is that the method satisfies the integral conservation of the quantities, such as mass, amount of linear movement and energy for any group of control volumes and, consequently, for the entire calculation domain [2].
The numerical solution methods of the system of algebraic equations given by Eq. (1) can be classified as direct or iterative [4]. Since the generated system of algebraic equation is sparse and large-scale, direct methods are not suitable due to their high computational cost. For this reason, we chose to work with iterative methods (solvers). Solvers are efficient in relaxing highfrequency (oscillatory) components during the first few iterations in refined grids. However, after some iterations, the process becomes slow, what means there is a predominance of low-frequency (smooth) modes [5,6].
Multigrid is an efficient method to accelerate solvers' convergence rate [5][6][7]. This method sweeps, during the iterative process, several grids of different refinement levels. According to Trottenberg et al. [6], smooth modes of the error become more oscillatory in coarser grids. Because of this, the different components of the error are efficiently smoothed, what accelerates the convergence of the iterative process used to solve the system given by Eq. (1).
Several studies concern the use of the multigrid method (MG) in problems involving a system of curvilinear coordinates. Smith et al. [8] presented a MG algorithm in order to accelerate a three-dimensional Navier-Stokes equations on general curvilinear grids. The authors calculated the metrics of the coarser grid while restricting the metrics of the adjacent finer grid. Moreover, they presented the convergence rate, computational times and performance indexes of the MG for several Reynolds numbers in highly curved ducts, with Full Approximation Scheme (FAS) combined with Full Multigrid (FMG).
Oosterlee et al. [9] found benchmark solutions for Navier-Stokes equations for incompressible flows, in steady state, in L-shaped domain, with systems of general curvilinear coordinates. The discretized system was solved using the MG with FAS, F-cycle, number of relaxations performed at the restriction and prolongation step (pre-and post-smoothing) equals to 1 and maximum number of levels.
Trottenberg et al. [6] stated that the choice of the components of the MG, such as solver, number of relaxations and type of cycle as well as restriction and prolongation operators, can strongly influence the convergence rate of the algorithm. According to Ferziger et al. [3], in the general context of the MG, many parameters can be chosen more or less arbitrarily, and the convergence rate and CPU time depend on these choices.
Li et al. [10] studied Navier-Stokes equations for two-dimensional flow of incompressible fluid as well as heat transfer in parallelogrammic, trapezoidal and sine-shaped cavities. The authors used FVM to discretize equations with collocated grid arrangement of the variables to solve physical problems in grids with 1024 × 1024 volumes. The system of algebraic equations was solved using the MG with FAS, V-cycle, and Strongly Implicit Procedure (SIP) solver as well as the SIMPLE method to calculate the pressure-velocity coupling. They presented numerical results for the problem of lid-driven cavity in parallelogrammic cavity and Reynolds numbers of 100 and 5000, more accurate than the solutions presented by [11]. Moreover, they presented benchmark solutions for natural and mixed-convection problems.
Two-dimensional Burgers' equations are regarded as one of the important models in the study of fluid flows. These equations have many applications, which go from cosmology to traffic modeling [12]. Ferm et al. [13], Zhang et al. [14], Neveu et al. [15] and Santiago et al. [16] studied the problem of Burgers' equations with the aim of improving the convergence of the MG. All of the authors used orthogonal structured grids. Among the studies found on MG applied to Burgers' equations, only a few aim at improving the convergence of the method by means of reducing the CPU time, and even less are systematic studies of the parameters of the MG.
In this paper, our goal is to systematically analyze some parameters of the geometric MG aimed at developing an efficient code for the problem of two-dimensional Burgers' equations, using systems of general curvilinear coordinates. In order to generate the non-orthogonal structured curvilinear grids, algebraic and differential methods were used [17]. Equations were discretized by FVM with Central Difference Scheme (CDS). The system of algebraic equations was solved with the solvers Lex-GS, η-line-GS, MILU and MSI. The metrics of the transformation and of the residual were restricted in order to solve the equations with curvilinear grids [8].
This paper is organized as follow: Section 2 presents the mathematical model in Cartesian and curvilinear coordinates; the generation of the coordinate system is shown in Section 3; Section 4 details the numerical model; Section 5 depicts the resolution methods for the system of linear equations; the MG is presented in Section 6; Section 7 presents the calculation of transformation metrics; Section 8 shows the results and discussion; and, finally, the conclusion is presented in Section 9.

Mathematical Model
Considering constant properties, steady state and Cartesian coordinates, the two-dimensional Burgers' equations are written dimensionless as and ∂uv ∂x where p stands as the static pressure, u and v are the velocity components on the coordinate directions x and y, respectively; and Re is the Reynolds number.
The source term S(x, y, Re) and the field pressure p are, respectively, In this study, we consider Re = 1. The analytical solutions, obtained by manufactured solution [18], are given by the expressions u(x, y) = f (x) g (y) and v(x, y) = −f (x)g(y). Thus and v (x, y) = −5x 4 + 10x 3 − 105 16 The computational domain of this problem, depicted in Fig. 1, is defined by In general, Eqs. (2) and (3) can be written as where φ (being u or v), p u = − ∂p ∂x and p v = − ∂p ∂y . In order to transform Eq. (6), written on the physical domain (x,y), to the computational domain (ξ ,η), the transformation of coordinates that is needed for two-dimensional problems is given by The metrics of the transformation are given by and the jacobian is given by Following the chain rule, the derivatives of the advective and diffusive flow of the generic variable ϕ and the terms of pressure, which appear in Eq. (6), result in ξ + y 2 ξ are the components of the metric tensor in two dimensions. After defining the following variables U = uy η − vx η and V = vx ξ − uy ξ , the advective flow can be rewritten as where U and V are the contravariant components of the velocity vector.
By grouping the transformed terms, we obtain Eq. (7) is the general equation transformed for a scalar φ.

Generation of the Coordinate System
The methods for grid generation are basically divided into algebraic and differential methods [17,19]. In this paper, we used an algebraic method that utilizes Lagrange interpolation and a differential method that uses the solution of elliptic equations for grid generation.

Numerical Model
By integrating Eq. (7) over the control volume P in the transformed plane, as seen in Fig. 3a, results in where L p φ J P is the approximation of the integral in ξ and η of p ϕ in the volume P, and w, e, s and n are the west, east, south and north faces of the volume P, respectively.
where λ e = 1 2 sign (U e ), λ w = 1 2 sign (U w ), λ n = 1 2 sign (V n ) and λ s = 1 2 sign (V s ), where sign stand as the sign function. This expression can be rewritten as (9) where NB stands as the 8 neighboring volumes of P.
In order to transform the conservation equations to the domain (ξ , η), it is necessary to know the transformation metrics, that is, the variables x η e P , y η e P , x ξ n P and y ξ n P , for the two-dimensional case. To make the notation easier, the face is used as superscript and the reference point as subscript. Fig. 3b shows a basic volume P, with points A, B, C and D in the respective vertex.
With these variables, it is possible to calculate the jacobian of the transformation. The metrics are numerically calculated by Therefore, the jacobian is calculated as where c is the center of the control volume and ψ ξ

Iterative Methods and Multigrid Details
A solution for Eq. (9) can be obtained by iterative method. In this section we will briefly describe those used in this article.

Gauss-Seidel Method and Its Variations
If the diagonal entries of A are non-zeros, the unknown corresponding value can be isolated in each equation, resulting in the lexicographic Gauss-Seidel method (Lex-GS) [4,20] and its variations.
It can be classified as a point-wise or block solver. In point-wise methods, each variable is updated individually. Fig. 4a depicts the lexicographic order, which was used in this study. By combining the Gauss-Seidel method with this order, we have the Lex-GS.
where the superscript m represents the mth iteration, the subscript the position on the grid, i = W , S, SW, SE and j = E, N, NW, NE. Notice that the nine closest neighbors were required in order to approximate ϕ p .
On block methods (columns or lines, for instance), each block is updated at once. Fig. 4b shows grid points ordered in lines in the direction η. By using the Gauss-Seidel method with this order, we have the η-line Gauss-Seidel (η-line-GS), which was used in this study.

Modified Incomplete LU Decomposition
Incomplete LU decomposition (ILU) of the A [20,21] consists in decomposing a matrix A in an incomplete manner. Such decomposition is of the A = LU − R form, where, in this subsection, L and U represent lower and upper triangular matrices and R represents the residue or decomposition error.
By doing the nine-point ILU decomposition of the coefficient matrix given by Eq. (9), we obtain L and U matrices with the same sparsity of A [20]. This case is called ILU(0).  This method is named Modified ILU decomposition (MILU) [6,20,22].
Given the MILU decomposition, the system given by Eq. Thus, the solution in the iteration m + 1 is given by φ (m+1) = φ (m) + e (m) .

Modified Strongly Implicit
In the MSI method [16,23], a LU decomposition is proposed, such that, M = LU, where the matrices L and U are lower and upper triangular, respectively, and the main diagonal of U is unitary. By performing the LU decomposition, in M four diagonals are filled. These non-null additional diagonals, are denoted by φ 1 i,j , φ 2 i,j , φ 3 i,j and φ 4 i,j . Therefore, the matrix M can be written as M = A + N, where A is the coefficient matrix of the equation system given by Eq. (9) and N consists only of diagonals φ 1 i,j , φ 2 i,j , φ 3 i,j and φ 4 i,j . A parameter σ was employed to partly cancel the influence of the additional terms that appear in M. By decomposing matrix A, the system is solved using the same iterative process shown in the Subsection 5.2.

Multigrid Details
Basic iterative methods are efficient in relaxing high-frequency components in refined grids during the first few iterations. However, after some iterations, the process becomes slow, what signals the predominance of low-frequency modes [5]. At this moment, it is recommended to transfer the information to the immediately coarser grid, where smooth error modes become more oscillatory and relaxation will be more efficient [6,22]. This happens because when the multigrid method (MG) is associated with an iterative method, it relaxes the error and corrects the solution in different grid sizes.
The MG sweeps a group of grids with different spacings. By means of a solver, iterations are performed at each level of the grid until the specified convergence criterion is reached for the system of equations of the finest grid. The sequence through which the grids are swept is denominated cycle. The coarsening ratio (r) of the grids is defined as r = H/h, where h is the size of the volumes of the finer grid and H represents the size of the volumes of the immediately coarser grid. The number of grids employed is called the number of levels (L). In case the highest possible number of grid levels is used, it will be denoted as L max .
The multigrid method can be implemented with Correction Scheme (CS), more suitable for linear problems, or Full Approximation Scheme (FAS), more indicated for non-linear problems [5]. In the FAS, the residual and the approximation of the solution are transferred to the coarser grids [6].
Operators that transfer information from a fine grid, h , to an immediately coarser grid, 2h , are called restriction operators (I 2h h ), the reverse are called prolongation operators (I h 2h ). The number of iteration used in the restriction and prolongation steps are called pre-and post-smoothing (ν 1 and ν 2 , respectively).

Calculation of Coordinates Transformation Metrics in Coarser Grids
It is also necessary to restrict the transformation metrics for the solution of the equations in curvilinear grids. In this paper, the convergence of the MG was achieved by calculating the coordinate transformation metrics as ψ η  The jacobian of the transformation on the coarse grid is calculated as where c is the center of the control volume, ψ ξ

Implementation Data
The algorithm was implemented in double precision FORTRAN 95 Intel 11.1 compiler. Simulations were performed in a 3.4 GHz Intel(R) Core(TM) i7-3770 microcomputer, 16 GB RAM, 64-bits Windows XP.
Geometric MG with FAS and V-cycle was used to solve the system of algebraic equations represented by Eq. (9) for the domain depicted in Fig. 1. The restriction operator of the approximations of the solution was obtained by arithmetic mean of the property values of the four volumes of the fine grid [6]. The restriction of the residual was done by adding the relative residuals of the control volumes of the fine grid that correspond to those of the coarser grid [24]. The correction prolongation was done by bilinear interpolation [5,6]. In this study, standard grid coarsening ratio was used (r = 2). Lex-GS, η-line-GS, MSI and MILU were used as solvers. The number of iterations in the pre-and post smoothing is the same, that is, ν = ν 1 = ν 2 . The number of unknowns is given by N = N ξ N η , where N ξ and N η stand for the number of volumes in the directions ξ and η, respectively. The stop criterion used is the l 1 norm of the residual in the current iteration r (m) , non-dimensionalized by the residual in the initial estimate, r (0) , that is ||r (m) || 1 /||r (0) || 1 ≤ tol , where tol = 10 −11 is adopted.

Average Convergence Factor
The empirical average convergence factor, which approximates the asymptotic convergence factor, is computed based on the residual. Such factor, as described in Trottenberg et al. [6], is given by where m represents the number of iterates or multigrid cycles performed.
In order to reduce the effect caused by the discard of elements during the ILU decomposition, a study of ρ m for different values of σ , with 0 ≤ σ ≤ 0.95 and −0.25 ≤ σ ≤ 1, for the MSI and MILU solvers, respectively, was carried out. For this, we considered L = L max and the number of inner iterations ν 1 = ν 2 = 3, for both solvers.
Figs. 7a and 7b depict the MSI method for the results of ρ m vs. σ for the variables u and v and grids generated by Lagrange interpolation and elliptic equations. Figs. 7c and 7d depict the same to MILU. Notice that the smallest values of ρ m , regardless of the size of the problem, happen when σ ≈ 0.9 and σ ≈ −0.2 for MSI and MILU solvers, respectively. Therefore, from now on, these values will be employed for σ in the formulations of MSI and MILU. Fig. 8 shows the empirical average convergence factor for all solvers assessed. The presented results are for different grid sizes generated by Lagrange interpolation and elliptic equations. As shown in Fig. 8, the η-line-GS method present the better empirical average convergence factor when compared with the Lex-GS method. This happened as the equations involve a strong coupling of the unknowns in the direction η.
As Trottenberg et al. [6] stated, stretched volumes in the physical domain cause anisotropy in the equations. In anisotropic problems, the convergence of basic methods, such as pointwise solvers, decreases [9,25,26]. The MILU and MSI solvers present the best empirical average convergence factors, that is, present ρ m ∼0 and ρ m 1, which is a desirable property. For instance, in the N = 4096 × 4096 grid, ρ m ≈ 0.03. Based on these results and aimed to develop an efficient algorithm to solve the problem in question, we will use the MILU and MSI solvers in our algorithm hereinafter as they have the best average convergence factors.

Inner Iterations
The effect of the number of inner iterations, denoted by ν = ν 1 = ν 2 , was assessed for grids sizes N = 512 × 512, 1024 × 1024, 2048 × 2048 and 4096 × 4096, using the maximum number of levels (L = L max ) in every case. Fig. 9 illustrate the influence of the number of inner iterations in the sum of the CPU times (t CPU ) in determining the velocities u and v, using the MSI and MILU solvers. In each curve, the value of ν that results in the lowest CPU time is indicated by the symbol . Noticeably, the lowest CPU time, for the algorithm built with the MSI solver, was obtained with four iterations, except in the grid N = 4096 × 4096, generated by Lagrange interpolation, and grid N = 2048 × 2048, generated by elliptic equations. Calculating the weighted average in CPU time gains for ν = 2 and ν = 4, the number of inner iterations of the solver that obtained the best average performance was ν = 4, for both grids generators. Thus, hereinafter, for the algorithm built with the MSI solver, the ν value adopted will be ν = 4.
For the algorithm built using MILU as solver, the lowest CPU time was obtained with three iterations for any of the values of N used, for grids generated by means of Lagrange interpolation. For grids generated by means of elliptic equations, the lowest CPU time was obtained with two iterations, except for the finest grid. Evaluating the weighted average in the CPU time gains for ν = 2 and 3, the number of inner iterations of the solver that obtained the best average performance was ν = 3, is approximately 1% lower. For this reason, the value adopted in the following studie is ν = 3, for the algorithm built using MILU as solver for both grids generators. One must notice that when the value of ν increases or decreases in relation to the ν adopted, the CPU time increases. Depending on the value of ν, this increase can be significant. The analysis also showed that the increase in CPU time is the same, or very close, for both grid generators, as shown in Fig. 9. Santiago et al. [16] solved two-dimensional Burgers' equations by means of the finite difference method in Cartesian coordinates, on orthogonal grids, for the square cavity problem, and tested MSI as solver, which was also used in this work. The authors found ν optimum = 5 using the MG with FAS (the authors did not mention the value of σ used).

Number of Grid Levels
Another important parameter of MG is the number of grid levels. For instance, considering a problem with r = 2, the coarsening ratio used in this study, and N = 128 × 128 unknowns, the highest possible number of grid levels is L max = 7 grids, which are 2 × 2, 4 × 4, 8 × 8, 16 × 16, 32 × 32, 64 × 64 and 128 × 128.
For the study of influence of the number of grid levels of MG over CPU time, we consider the number of inner iterations suggested in the previous section. Fig. 10 depict the effect of the number of grid levels in the sum of the CPU time in order to determine the velocities u and v, with MSI and MILU as solvers. The symbol indicates the L that resulted on the lowest CPU time (L optimum ) in each curve. Fig. 10a (Lagrange interpolation) that the number of grid levels in which the lowest CPU time was obtained is L optimum = L max − (1, 2 or 3). For grids generated by employing elliptic equations (Fig. 10b), L optimum = L max − (1, 3 or 4). In Fig. 10c (Lagrange interpolation), it is possible to observe that the number of grid levels in which the lowest CPU time was obtained is L optimum = L max − (1, 2 or 3). For grids generated by employing elliptic equations, Fig. 10d, is L optimum = L max − (0 or 1). Moreover, Fig. 10 show that for both grid generators, the CPU time can significantly increase according to the number of levels used. This is the case for small number of grid levels, that is, number of grid levels close to L = 1, which is the case of the singlegrid method. However, for L ≡ L max , with MSI or MILU as solvers, the CPU time does not change significantly. It was found that the number of grid levels what present the best average performance (weighted average of the percentages in gains in CPU time) was L = L max − 1 for algorithms built using MSI as solver, for the four biggest grid sizes and for both grid generators used. For MILU solver, the best average performance was obtained for L = L max − 3 for grids generated using Lagrange interpolation and L = L max for grids generated by employing elliptic equations, considering the four biggest grid sizes.

Notice in
By means of the FVM, Rabi et al. [27] solved the conductive-convective problem with the multigrid method by employing orthogonal, structured grids for a problem with N = 64 × 64 volumes and recommend using at least L = 4 for good error relaxation. Kumar et al. [24] state that there is no gain with L higher than 4 for the lid-driven cavity problem using orthogonal structured grids, FVM and N = 513 × 513 volumes. For Burgers' equations using orthogonal structured grids and MG with FAS, Santiago et al. [16] reached L optimum = L max − (0 to 3) by employing FDM. The authors considered L optimum = L max , since for L optimum < L ≤ L max the CPU time was virtually the same.  Fig. 11 for the MSI and MILU solvers.
As shown in Fig. 11, the MG was extremely efficient. According to Ferziger et al. [3], the higher the number of unknowns N, the higher the advantage of MG over SG. Notice that the theoretical performance of the MG for Burgers' equations using the general curvilinear coordinates system was confirmed, p ≈ 1. Besides, p ≈ 2 for the SG, what was also expected [4]. Furthermore, the CPU time and p values were very close for both methods, with a slight advantage of the MSI.
Tab. 2 shows execution time, in seconds, of the SG and MG with MSI and MILU as solvers as well as the grid generators using Lagrange interpolation and elliptic equations. Based on the data presented in Tab. 2, considering MG with FAS, the CPU times for the problem with the MSI solver and grids generated using Lagrange interpolation, are smaller than those obtained with MILU, being approximately 19% faster in the most refined grid. For the problem with the MSI solver and grids generated using elliptic equations are also smaller than those obtained with MILU, being approximately 16% faster in the most refined grid. Comparing SG CPU times, in Tab. 2, MSI was also faster in all grids. Tab. 3 shows the discretization error in infinity norm of the SG and MG with MSI as solver as well as the grid generators using Lagrange interpolation and elliptic equations.
Considering the results of the Tab. 3, the solutions obtained for the problem with grids generated using elliptic equations are slightly more accurate than those using Lagrange interpolation. In this sense, the MG with MSI as solver and grids generated by employing Lagrange interpolation stood out as the best combination (lowest CPU time and discretization error) when compared to the other possible combinations studied.

Conclusions
In this paper, two-dimensional Burgers' equations in steady state with constant properties, using a system of curvilinear coordinates in an L-shaped domain was solved. The influence of some parameters of the geometric multigrid method as well as some solvers to solve this problem were assessed. In order to generate grids, Lagrange interpolation and systems of elliptic equations were used. Equations were discretized by FVM, generating systems of linear equations, which were solved by applying the Geometric Multigrid method and the GS-Lex, η-line-GS, MSI and MILU solvers. Convergence of the multigrid method was achieved by calculating the coordinate transformation metrics as presented in this paper. Based on the obtained results, it was mainly found that from the solvers employed, MSI and MILU are the most efficient for the problem assessed, with empirical average convergence factor close to 0.03; results show that the MSI solver, with ν = 4 and L = L max − 1, is faster than MILU for both grid generators, since it has the best complexity factors values. In the grid N = 4096 × 4096, the multigrid method with MSI as solver is approximately 19% e 16% faster than the MILU solver, for grids generated by Lagrange interpolation and elliptic equations, respectively. The algorithm built using the MG with MSI as solver and grids generated by employing Lagrange interpolation stood out as the best combination (lowest CPU time and discretization error) when compared to the other possible combinations studied.
It is also necessary to restrict the coordinates transformation metrics to solve this problem in curvilinear grids. Therefore, it shows how to calculate the coordinates transformation metrics in the coarser grids.
acknowledge the infrastructure of the Numerical Experimentation Laboratory (LENA) and the Graduate Program in Numerical Methods in Engineering of the Federal University of Parana (UFPR).

Funding Statement:
The author(s) received no specific funding for this study.

Conflicts of Interest:
The authors declares that they have no conflicts of interest to report regarding the present study.